timingAnalyzer.cpp

#include "TSystem.h"
#include "TApplication.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TGraph.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TProfile.h"
#include "TChain.h"
#include <TTreeReader.h>
#include <TTreeReaderValue.h>
#include <TTreeReaderArray.h>
//#include <ProcessorRootStruc.hpp>
#include <PaassRootStruct.hpp>
#include <cmath>
#include <tuple>
#include "position.h"

std::tuple<double,double,double> pcfdAnalyzer(vector<double> trace, double frac);
double regCfd(vector<double> &trace);
double derAnalyzer(vector<double> &trace);
pair <double, double> baselineCalc(vector<double> trace);
double traceAnalyzer(vector<double> trace);
double fittingAnalyzer(vector<double> trace);
double gFit(vector<double> trace);
vector<double> reducer(vector<double> trace);
double centroid(vector<double> trace);
TF1 *fit = new TF1("fit","(x>[1])*[0]*exp(([1]-x)/[2])*(1-exp(-1*pow(x-[1],[5])/[3]))+[4]",0,300);
TF1 *fpol1 = new TF1("fName1","pol1",0,300);
TF1 *fpol3 = new TF1("fName3","pol3",0,1);
TF1 *fpol2 = new TF1("fName2","pol2",0,1);
TF1 *gauss = new TF1("guass","gaus",0,1);

void timingAnalyzer(){
  TFile *_file0 = TFile::Open("ysoSingle_DD.root");
  TTree *GS = (TTree*)_file0->Get("PixTree");
  TTreeReader singe;
  cout<<"readermade"<<endl;
  singe.SetTree( GS );
  //TTreeReaderArray<processor_struct::ROOTDEV> rd = {singe, "root_dev_vec_"};  //gives vector of stucture 
  TTreeReaderArray<processor_struct::ROOTDEV> rd = {singe, "rootdev_vec_"};  //gives vector of stucture 

    TH1D *plot1=new TH1D("plot1","plot1",4000.,-100.,100.);
    TH2D *energy=new TH2D("energy","energy",5000.,0.,10000.,5000.,0.,10000.);
    TH2D *phase=new TH2D("phase","phase",1000.,80.,105.,1000.,80.,105.);

    //variables
    int eventNum = 0;

    double energyone=0.0, energytwo=0.0, energythree=0.0;
    vector<double> traceone, tracetwo, tracethree, xatrace, xbtrace, yatrace, ybtrace;
    vector<unsigned int> tone, ttwo, tthree, xat, xbt, yat, ybt;
    int ione=0, itwo=0, ithree=0, twocount=0;
    double hrTimeone=0.0, hrTimetwo=0.0, hrTimethree=0.0;
    double cfdTimeone=0.0, cfdTimetwo=0.0, cfdTimethree=0.0;
    double phaseone=0.0, phasetwo=0., tmaxone=0., tmaxtwo=0.;
    double timeone=0.0, timetwo=0.0, timethree=0.0;
    double qdcone=0.0, qdctwo=0.0;
    double tdiff=0.0, baselineone=0.0, baselinetwo=0.0, sigmaone=0., sigmatwo=0.;
    double ecaltwo=0.0, slopeone=0., slopetwo=0.;
    int maxvalone=0, maxvaltwo=0;
    double xa=0., xb=0., ya=0., yb=0.;
    double xaqdc=0., xbqdc=0., yaqdc=0., ybqdc=0.;
    double xatm=0., xbtm=0., yatm=0., ybtm=0.;
    double xposq=0., yposq=0.;
    double xpost=0., ypost=0.;

  TFile *newFile = new TFile("ysoTraces.root","RECREATE");
  TTree *timing = new TTree("timing","tree filled with traces and energies etc.");

  timing->Branch("eventNum", &eventNum);
  timing->Branch("energyone", &energyone);
  timing->Branch("energytwo", &energytwo);
  timing->Branch("energythree", &energythree);
  timing->Branch("energycal", &ecaltwo);
  timing->Branch("traceone", &traceone);
  timing->Branch("tracetwo", &tracetwo);
  timing->Branch("tracethree", &tracethree);
  timing->Branch("ione", &ione);
  timing->Branch("itwo", &itwo);
  timing->Branch("hrTimeone", &hrTimeone);
  timing->Branch("hrTimetwo", &hrTimetwo);
  timing->Branch("cfdTimeone", &cfdTimeone);
  timing->Branch("cfdTimetwo", &cfdTimetwo);
  timing->Branch("phaseone", &phaseone);
  timing->Branch("phasetwo", &phasetwo);
  timing->Branch("slopeone", & slopeone);
  timing->Branch("slopetwo", & slopetwo);
  timing->Branch("timeone", &timeone);
  timing->Branch("timetwo", &timetwo);
  timing->Branch("tmaxone", &tmaxone);
  timing->Branch("tmaxtwo", &tmaxtwo);
  timing->Branch("qdcone", &qdcone);
  timing->Branch("qdctwo", &qdctwo);
  timing->Branch("tdiff", &tdiff);
  timing->Branch("xposq", &xposq);
  timing->Branch("yposq", &yposq);
  timing->Branch("xpost", &xpost);
  timing->Branch("ypost", &ypost);
  timing->Branch("bone", &baselineone);
  timing->Branch("btwo", &baselinetwo);
  timing->Branch("sone", &sigmaone);
  timing->Branch("stwo", &sigmatwo);

  fit->SetParLimits(0,500,50000);
  fit->SetParLimits(2,5,70);
  fit->SetParLimits(3,40,200);
  //fit->SetParLimits(5,1.0,1.8); //power for fit
  fit->FixParameter(5,4.);

  std::vector<unsigned int> *trace;
  while(singe.Next()){

    if(eventNum%1000==0) cout<<"\r" << eventNum<< flush;
    tone.clear(); ttwo.clear(); tthree.clear();
    traceone.clear();
    tracetwo.clear(); tracethree.clear();
    xat.clear(); xbt.clear(); yat.clear(); ybt.clear();
    xatrace.clear(); xbtrace.clear();
    yatrace.clear(); ybtrace.clear();
    ione=0; itwo=0; ithree=0;
    hrTimeone=0.0; hrTimetwo=0.0; hrTimethree=0.0;
    energyone= -1000.0; energytwo= -1000.0; energythree= -1000.0;
    cfdTimeone=0.0; cfdTimetwo=0.0; cfdTimethree=0.;
    phaseone=0.0; phasetwo=0.0; slopeone=0.; slopetwo=0.;
    timeone=0.0; timetwo=0.0, timethree=0.0;
    tmaxone=0.0; tmaxtwo=0.0; tdiff=0.0;
    ecaltwo=0.0; baselineone=0.0; baselinetwo=0.0; sigmaone=0; sigmatwo=0;
    qdcone=0.0;qdctwo=0.0;
    maxvalone=0;maxvaltwo=0;
    xa=0; xb=0.; ya=0.; yb=0.;
    xposq=0.; yposq=0.;
    xpost=0.; ypost=0.;
    xaqdc=0.; xbqdc=0.; yaqdc=0.; ybqdc=0.;
    xatm=0.; xbtm=0.; yatm=0.; ybtm=0.;

    for(auto itC = rd.begin(); itC!=rd.end();itC++){
        trace = &(itC->trace);
        double energy = itC->energy;
        double time = itC->timeSansCfd;
        double highResTime = itC->highResTime;
        std::string dtype = itC->subtype.Data();
        std::string dgroup = itC->group.Data();
        int channel = itC->chanNum; 
        int module = itC->modNum;
        int maxValue = itC->maxVal;
        if(eventNum==0) cout<<trace->size()<<" "<<channel<<endl;

        if (channel==0 && module==1){ 
            ione++;
            tone = *trace;
            double hold =0.0;
            for(int i=0;i<tone.size();i++){
              hold = (double) tone[i];
              traceone.push_back(hold);
              if(traceone[tmaxone]<=traceone[i]) tmaxone=i;
            }
            //for(int i=0;i<traceone.size()-2;i++)traceone[i] = (traceone[i]+traceone[i+1])/2;
            energyone = energy;
            hrTimeone = highResTime;
            cfdTimeone = time;
            timeone = time;
            maxvalone = traceone[tmaxone];
            //phaseone = cdfAnalyzer(traceone);
            //cfdTimeone += 2*phaseone;
            qdcone = traceAnalyzer(traceone);
            auto baselineVals = baselineCalc(traceone);
            baselineone = baselineVals.first;
            sigmaone = baselineVals.second;
        }

        else if(channel==12 && module==0){ //external trig
            itwo++;
            ttwo = *trace;
            double hold =0;
            for(int i=0;i<ttwo.size();i++){
              hold = (double) ttwo[i];
              tracetwo.push_back(hold);
              if(tracetwo[tmaxtwo]<=tracetwo[i]) tmaxtwo=i;
            }
            //for(int i=0;i<tracetwo.size()-2;i++)tracetwo[i] = (tracetwo[i]+tracetwo[i+1])/2;
            energytwo = energy;
            hrTimetwo = highResTime;
            cfdTimetwo = time;
            timetwo = time;
            maxvaltwo=maxValue;
            //phasetwo = cdfAnalyzer(tracetwo);
            //cfdTimetwo += 2*phasetwo;
            qdctwo = traceAnalyzer(tracetwo);
            auto baselineVals = baselineCalc(tracetwo);
            baselinetwo = baselineVals.first;
            sigmatwo = baselineVals.second;
        }
        else if(channel==13 && module==0){ //hagrid
            ithree++;
            tthree = *trace;
            double hold =0;
            for(int i=0;i<tthree.size();i++){
              hold = (double) tthree[i];
              tracethree.push_back(hold);
            }
            energythree = energy;
            //qdcone = traceAnalyzer(tracethree);
            hrTimethree = highResTime;
            //cfdTimethree = time;
        }

        else if(channel==0 && module==1){
          xat = *trace;
          double hold =0;
          for(int i=0;i<xat.size();i++){
              hold = (double) xat[i];
              if(hold>xatm) xatm = hold;
              xatrace.push_back(hold);
            }
          xa = energy;
          xaqdc = traceAnalyzer(xatrace);
        }
        else if(channel==1 && module==1){ 
          xbt = *trace;
          double hold =0;
          for(int i=0;i<xbt.size();i++){
              hold = (double) xbt[i];
              if(hold>xbtm) xbtm = hold;
              xbtrace.push_back(hold);
            }
          xb = energy;
          xbqdc = traceAnalyzer(xbtrace);
        }
        else if(channel==2 && module==1){ 
          yat = *trace;
          double hold =0;
          for(int i=0;i<yat.size();i++){
              hold = (double) yat[i];
              if(hold>yatm) yatm = hold;
              yatrace.push_back(hold);
            }
          ya = energy;
          yaqdc = traceAnalyzer(yatrace);
        }
        else if(channel==3 && module==1){ 
          ybt = *trace;
          double hold =0;
          for(int i=0;i<ybt.size();i++){
              hold = (double) ybt[i];
              if(hold>ybtm) ybtm = hold;
              ybtrace.push_back(hold);
            }
          yb = energy;
          ybqdc = traceAnalyzer(ybtrace);
        }
    }
    if((ione>0 && itwo>0) || maxvalone>4000){
      std::tuple<double,double,double> result = pcfdAnalyzer(traceone,0.1);
      //phaseone = 4*fittingAnalyzer(traceone);
      phaseone = 4*(std::get<0>(result));
      tmaxone = std::get<1>(result);
      slopeone = std::get<2>(result);
      cfdTimeone += phaseone;
      result = pcfdAnalyzer(tracetwo,0.1);
      phasetwo = 4*(std::get<0>(result));
      tmaxtwo = std::get<1>(result);
      slopetwo = std::get<2>(result);
      cfdTimetwo += phasetwo;
      tdiff = cfdTimeone-phaseone-cfdTimetwo+phaseone;
      ecaltwo = 0.09369787313*energytwo-1.639424969;
      pair<double,double> posq = position(xaqdc,xbqdc,yaqdc,ybqdc);
      pair<double,double> post = position(xa,xb,ya,yb);
      xposq = posq.first;
      yposq = posq.second;
      xpost = post.first;
      ypost = post.second;
      //plot1->Fill(tdiff);
      //energy->Fill(energyone,energytwo);
      //phase->Fill(phaseone,phasetwo);
      timing->Fill();
    }
    eventNum++;
  }
  timing->Write();
  //plot1->Write();
  //phase->Write();
  //energy->Write();
  cout<<endl;
}

std::tuple<double,double,double> pcfdAnalyzer(vector<double> trace, double frac){
  int maxPos=0;
  vector<double>::iterator it;
  TGraph *fTraces = new TGraph();
  double fmax=0.0;
  double pmax=0.0;
  double baseline=834.;
  double thresh=0.0;
  double phase=0.0;
  std::pair <double,double> points(0.,0.);

  if(trace.size()>10){
    it = max_element(trace.begin(),trace.end());
    maxPos = distance(trace.begin(),it);
  }
  else{
    return make_tuple(-9999.,-9999.,-9999.);
  }
  for (int iv = 0;iv<trace.size();iv++){
    fTraces->SetPoint(iv, iv, trace.at(iv));
  }
  
  std::pair <double, double> range((maxPos-2),(maxPos+1));
  fpol3->SetRange(range.first,range.second);
  fTraces->Fit(fpol3,"RNQSW");
  if(!(fpol3->IsValid()))return make_tuple(-5555.,-5555.,-5555.);
  fmax = fpol3->GetMaximum(range.first,range.second);
  pmax = trace.at(maxPos);
  //if(frac>0.11)
  baseline = baselineCalc(trace).first;
  if(pmax<3*baselineCalc(trace).second)cout<<"peak too small "<<baselineCalc(trace).second<<endl;
  //if(maxPos>125)return -7777.;
  thresh = (fmax-baseline)*frac+baseline; //here frac is the Constant Fraction
  for(int i=maxPos;i>0;i--){
    if((trace.at(i)>=thresh) && (trace.at(i-1)<thresh)){
      points.first = i-1;
      points.second = i; //setting points of trace around threshold 
    }
  }
  fpol1->SetRange(points.first,points.second);
  fTraces->Fit(fpol1,"RNQSW");
  phase = fpol1->GetX(thresh,points.first,points.second);
  return make_tuple(phase,fmax-baseline,fpol1->GetParameter(1));
}

double regCfd(vector<double> &trace){
  double phase=0.0;
  vector<double> invTrace(trace.size(),0.0);
  double baseline=0.0;
  baseline = baselineCalc(trace).first;
  int highbound=0, lowbound=1000;
  int maxPos=0;
  vector<double>::iterator it;
  vector<double>::iterator itMin;
  vector<double>::iterator itMax;
  TGraph *fTraces = new TGraph();

  if(trace.size()>10){
    it = max_element(trace.begin(),trace.end());
    maxPos = distance(trace.begin(),it);
  }
  else
    return -9999.; // This is about the range of the waveform

  double totMin=0;
  double totMax=0;
  double zeroPos=0;

  //if(trace[maxPos]-baseline<100)
    //return 600;

  for(int i=0;i<trace.size();i++){
    trace[i] -= baseline;
    invTrace[i]=(-0.7*trace[i+3]);
    trace[i]+=invTrace[i];
  }
  for(int i=maxPos;i>maxPos-3;i--){
    if(trace[i]<0){
      zeroPos=i;
      break;
    }
  }

  for (int iv = 0;iv<trace.size();iv++){
    fTraces->SetPoint(iv, iv, trace.at(iv));
  }
  
  double m=0.0, b=0.0;
  m = (trace.at(zeroPos+1)-trace.at(zeroPos))/1;
  b = trace.at(zeroPos+1)-m*(zeroPos+1);
  phase = -1*b/m; 

  if(phase<0.0) return -1000;

  /*fpol1->SetRange(zeroPos,zeroPos+1);
  fTraces->Fit(fpol1,"RNQSW");
  phase = fpol1->GetX(0.0,zeroPos,zeroPos+1);*/
  return phase;
  
}

double derAnalyzer(vector<double> &trace){
  int maxPos=0;
  vector<double>::iterator it;
  TGraph *fTraces = new TGraph();
  TGraph *fTraces1 = new TGraph();
  //TF1 *fpol2 = new TF1("fName2","pol2",0,1);
  double fmax=0.0;
  double phase=0.0;
  double par[2]={0.0,0.0};
  vector<double> der;
  double hold=0.0;

  if(trace.size()<=10){
    return -9999.;
  }

  for (int iv = 0;iv<trace.size();iv++){
    fTraces1->SetPoint(iv, iv, trace.at(iv));
  }
  der.push_back(0.0);
  for(int i=1;i<trace.size()-1;i++){
    /*fpol1->SetRange(i-1,i+1);
    fTraces1->Fit(fpol1,"RNQSW");
    fpol1->GetParameters(&par[0]);
    der.push_back(par[1]);*/
    hold = (trace[i+1]-trace[i-1])/2.0;
    der.push_back(hold);
  }
  der.push_back(0.0);
  /*for(int i=0;i<trace.size();i++){
    //if(der[i]>der[maxPos]) maxPos=i;
    trace[i] = der[i];
  }*/
  it = max_element(der.begin(),der.end());
  maxPos = distance(der.begin(),it);
  for (int iv = 0;iv<der.size();iv++){
    fTraces->SetPoint(iv, iv, der.at(iv));
  }

  pair <double, double> range((maxPos-1),(maxPos+1));
  fpol2->SetRange(range.first,range.second);
  fTraces->Fit(fpol2,"RNQSW");
  fmax = fpol2->GetMaximum(range.first,range.second);
  phase = fpol2->GetX(fmax,range.first,range.second);
  return phase;

}

pair <double, double> baselineCalc(vector<double> trace){
    //calculating the baseline
    const int numBins = 20;
    double baseSum = 0.;
    for(int j=0;j<numBins;j++){
        baseSum += trace.at(j);
    }
    double  baseline = baseSum/(double(numBins));
    
    //calculating the standard dev
    double stddev = 0.0;
    for(int j=0;j<numBins;j++){
        stddev += pow(trace.at(j)-baseline,2);
    }
    stddev = sqrt(stddev/(double(numBins)));
    
    return std::make_pair (baseline, stddev);
}

double traceAnalyzer(vector<double> trace){
    double qdc=0.0;
    unsigned int highbound, lowbound;
    double baseline = baselineCalc(trace).first;
    int maxPos=0;
    vector<double>::iterator it;

    if(trace.size()>10){
    it = max_element(trace.begin(),trace.end());
    maxPos = distance(trace.begin(),it);
    }
    else if(maxPos<180){
      return -4444.;
    }
    else{
      return -9999.;
    }

    highbound = maxPos+40;
    lowbound = maxPos-20;

    for(int i=lowbound;i<highbound;i++)
      qdc += 0.5 * double(abs(trace[i-1] - baseline) + abs(trace[i] - baseline));

    return qdc;
}

double fittingAnalyzer(vector<double> trace){
  int maxLoc=0,dMax=0;
  double phase=0.0;
  if(trace.size()<10) return -9999.;
  double baseline=baselineCalc(trace).first;
  //TGraph *fTraces = new TGraph();
  TH1D *hist = new TH1D("hist","trace",250,0.,250.);
  double fmax=0.0;
  double der=0.0;
  
  double thresh=0.4;
  std::pair <double,double> points(0,0);


  //need to determine fitting range..want to start around beginning of rise
  for(int i=0;i<trace.size();i++){
    if(trace.at(i)>trace.at(maxLoc)) maxLoc=i;
    //fTraces->SetPoint(i, i, trace.at(i));
    if(trace[i+1]-trace[i]>der && i<=maxLoc && i<trace.size()-1){
      dMax=i; 
      der = trace[dMax+1]-trace[dMax];
    }
    for(int k=0;k<(int)trace.at(i);k++) hist->Fill(double(i));
  }

  /*for(int i=0;i<trace.size()-1;i++){
    //this is taking the derivative to determine the max risetime, which should be on the initial kink up
    if(trace[i+1]-trace[i]>der && i<=maxLoc)dMax=i; 
    der = trace[dMax+1]-trace[dMax];
  }*/
  for(int j=0;j<trace.size()-1;j++){
    double hold=0., temp=0.;
    temp = trace.at(j+1)-trace.at(j);
    if(temp>0) hold = temp;
    else  hold = -0.2*temp;
    //fTraces->SetPoint(j, 2*j, hold);
    double weight = hold/der+0.001;
    //fTraces->SetPoint(j, 2*j,weight);
    for(int k=0;k<(int)trace.at(j);k++) hist->Fill(double(j));
    hist->SetBinError(j,1/weight);
  }
  
  //general ranges for parameters, to help the fit actually work.
  //fit->SetParameter(1,dMax);
  fit->SetParLimits(0,500,50000);
  fit->SetParLimits(1,dMax-12,dMax);
  fit->SetParLimits(2,1,100);
  fit->SetParLimits(3,1,2000);
  fit->FixParameter(4,baseline);
  fit->FixParameter(5,3);

  std::pair <double, double> range(dMax-20,maxLoc+10);
  fit->SetRange(range.first,range.second);
  hist->Fit(fit,"RNSQW");
  /*fmax = fit->GetMaximum(range.first,range.second);
  //phase = fit->GetX(fmax,range.first,range.second);
  //phase = fit->GetParameter(1);
  phase = fit->GetParameter(1);
  //return phase;
  thresh = (fmax-baseline)*0.5+baseline; //here 0.5 is the constant fraction
  for(int i=maxLoc;i>0;i--){
    if((trace.at(i)>=thresh) && (trace.at(i-1)<thresh)){
      points.first = i-1;
      points.second = i;
    }
  }
  fpol1->SetRange(points.first,points.second);
  fTraces->Fit(fpol1,"RNSQW");
  phase = fpol1->GetX(thresh,points.first,points.second);
*/
  delete hist;
  return fit->GetParameter(1);
  //return phase;
}

double gFit(vector<double> trace){
  int maxPos=0;
  vector<double>::iterator it;
  TGraph *fTraces = new TGraph();
  double fmax=0.0;
  double pmax=0.0;
  double baseline=0.0;
  double phase=0.0;
  std::pair <double,double> points(0,0);

  if(trace.size()>10){
    it = max_element(trace.begin(),trace.end());
    maxPos = distance(trace.begin(),it);
  }
  else{
    return -9999.;
  }
  for (int iv = 0;iv<trace.size();iv++){
    fTraces->SetPoint(iv, iv, trace.at(iv));
  }
  
  std::pair <double, double> range((maxPos-2),(maxPos+2));
  gauss->SetRange(range.first,range.second);
  fTraces->Fit(gauss,"RNQSW");
  if(!(gauss->IsValid()))return -5555.;
  fmax = gauss->GetMaximum(range.first,range.second);
  phase = gauss->GetX(fmax,range.first,range.second);
  return phase;
}

vector<double> reducer(vector<double> trace){
  int n;
  n = trace.size();
  vector<double> t2;
  for(int i=0;i<n;i++){
    if(i%2==0) t2.push_back(0.5*(trace[i]+trace[i+1]));
  }
  return t2;
}

double centroid(vector<double> trace){

  double cent=0.0, ampSum=0.0, time=0.0;
  double baseline = baselineCalc(trace).first;
  double stddev = baselineCalc(trace).second;
  int maxLoc=0;
  pair <int, int> wave = std::make_pair(0,0);

  for(int i=0;i<trace.size();i++){
    trace[i]-=baseline;
    if(trace[i]>trace[maxLoc]) maxLoc=i;
  }
  for(int i=maxLoc-4;i<trace.size();i++){
    if(trace[i]*i>cent)cent=trace[i]*i;
    ampSum+=trace[i];
  }

  return cent/ampSum;
}