avgefficiency.C

Go to the documentation of this file.

// Macro for making plots of average cell efficiencies by view
// Runs on the output histogram file from HitEfficiency_module.cc,
// though you will want to either TChain or hadd a bunch of these
// files together to get any sort of reasonable statistics

#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"
#include "TLegend.h"
#include "TProfile.h"

#include <cmath>
#include <string>
#include <iostream>


void avgefficiency()
{

  double nmincell = 500; // mininum number of entries in denominator of cell efficiency for a cell to be "high" stat
  
  // Get the ana histo file and grab the relevant histograms
  TFile* tfile = new TFile("nova/ana/users/raddatz/raddatz/LowThreshold/haddlowthres.root","READ");
  
  // numerator of efficiency vs w (lencut histos are for path lengths 6-7 cm)
  TH2F* hitsVwcell = (TH2F*)tfile->Get("hiteff/hitsVwcelllencut");
  // denominator of efficiency vs w (lencut histos are for path lengths 6-7 cm)
  TH2F* cellsVwcell = (TH2F*)tfile->Get("hiteff/cellsVwcelllencut");
  // numerator of efficiency vs path length vs w for x cells
  TH2F* hitsVwVlenx = (TH2F*)tfile->Get("hiteff/hitsVwVlenx");
  // denominator of efficiency vs path length vs w for x cells
  TH2F* cellsVwVlenx = (TH2F*)tfile->Get("hiteff/cellsVwVlenx");
  // numerator of efficiency vs path length vs w for y cells
  TH2F* hitsVwVleny = (TH2F*)tfile->Get("hiteff/hitsVwVleny");
  // denominator of efficiency vs path length vs w for y cells
  TH2F* cellsVwVleny = (TH2F*)tfile->Get("hiteff/cellsVwVleny");
  
  // Define a bunch of histogram binning
  // First grab predefined binning from histograms
  int ncellbins = hitsVwcell->GetNbinsX();
  double cellmax = hitsVwcell->GetXaxis()->GetXmax();
  int nwbins = hitsVwcell->GetNbinsY();
  double wmin = hitsVwcell->GetYaxis()->GetXmin();
  double wmax = hitsVwcell->GetYaxis()->GetXmax();
  int npathbins = hitsVwVlenx->GetNbinsY();
  double pathmin = hitsVwVlenx->GetYaxis()->GetXmin();
  double pathmax = hitsVwVlenx->GetYaxis()->GetXmax();
  int neffbins = 101;
  double effmin = -0.05;
  double effmax = 1.05;
  double plmin = -0.5;
  int cellmaxplanecell = ((int)cellmax-1)%384;
  double plmax = ((int)cellmax-1-cellmaxplanecell)/384 + 0.5;
  int nplbins = (int)(plmax-plmin);
  int nclbins = 385;
  double clmin = -0.5;
  double clmax = 384.5;

  
  // Define a bunch of histograms 
  // efficiency as a function of path length and w histograms as well as ratio between views  
  TH2F* effVwVlenx = new TH2F("effVwVlenx",";W (cm); Path Length (cm)",nwbins,wmin,wmax,npathbins,pathmin,pathmax);
  TH2F* effVwVleny = new TH2F("effVwVleny",";W (cm); Path Length (cm)",nwbins,wmin,wmax,npathbins,pathmin,pathmax);
  TH2F* effVwVlenratio = new TH2F("effVwVlenratio",";W (cm);Path Length (cm)",
                                   effVwVlenx->GetNbinsX(),effVwVlenx->GetXaxis()->GetXmin(),effVwVlenx->GetXaxis()->GetXmax(),
                                   effVwVlenx->GetNbinsY(),effVwVlenx->GetYaxis()->GetXmin(),effVwVlenx->GetYaxis()->GetXmax());
   
  // efficiency vs. w averaged over cells for each view
  TH1F* xeffVw = new TH1F("xeffVw",";W (cm);Efficiency",nwbins,wmin,wmax);
  TH1F* yeffVw = new TH1F("yeffVw",";W (cm);Efficiency",nwbins,wmin,wmax);
  
  // efficiency vs w averaged over cells for each view using only cells with high stats
  TH1F* xeffVwcut = new TH1F("xeffVwcut",";W (cm);Efficiency",nwbins,wmin,wmax);
  TH1F* yeffVwcut = new TH1F("yeffVwcut",";W (cm);Efficiency",nwbins,wmin,wmax);
   
  // average efficiency Vs Cell Vs Plane for the lowest third, mid third, and highest third of the cell for x cells
  TH2F* effslowCellVPlanex = new TH2F("effslowCellVPlanex",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
  TH2F* effsmidCellVPlanex = new TH2F("effsmidCellVPlanex",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
  TH2F* effshighCellVPlanex = new TH2F("effshighCellVPlanex",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
  
  // averaged efficiency of the lowest third, mid third, and highest third of the cell for x cells
  TH1F* effslowx = new TH1F("effslowx",";Efficiency;NCells",neffbins,effmin,effmax);
  TH1F* effsmidx = new TH1F("effsmidx",";Efficiency;NCells",neffbins,effmin,effmax);
  TH1F* effshighx = new TH1F("effshighx",";Efficiency;NCells",neffbins,effmin,effmax);
  
  // average efficiency Vs Cell Vs Plane for the lowest third, mid third, and highest third of the cell for y cells
  TH2F* effslowCellVPlaney = new TH2F("effslowCellVPlaney",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
  TH2F* effsmidCellVPlaney = new TH2F("effsmidCellVPlaney",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
  TH2F* effshighCellVPlaney = new TH2F("effshighCellVPlaney",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
  
  // averaged efficiency of the lowest third, mid third, and highest third of the cell for y cells
  TH1F* effslowy = new TH1F("effslowy",";Efficiency;NCells",neffbins,effmin,effmax);
  TH1F* effsmidy = new TH1F("effsmidy",";Efficiency;NCells",neffbins,effmin,effmax);
  TH1F* effshighy = new TH1F("effshighy",";Efficiency;NCells",neffbins,effmin,effmax);
  
  // helper histogram to get the average efficiency of a cell over varius w ranges
  TH1F* eff = new TH1F("eff","",nwbins,wmin,wmax);
   
  int lowwbin = (int)((double)nwbins/3.0);
  int highwbin = (int)(2.0*(double)nwbins/3.0);
  // helper histograms for getting the average and variation in average cell efficiencies by view
  TH2F* effVwx = new TH2F("effVwx",";W (cm);Efficiency",nwbins,wmin,wmax,neffbins,effmin,effmax);
  TH2F* effVwy = new TH2F("effVwy",";W (cm);Efficiency",nwbins,wmin,wmax,neffbins,effmin,effmax);
  TH2F* effVwxcut = new TH2F("effVwxcut",";W (cm);Efficiency",nwbins,wmin,wmax,neffbins,effmin,effmax);
  TH2F* effVwycut = new TH2F("effVwycut",";W (cm);Efficiency",nwbins,wmin,wmax,neffbins,effmin,effmax);


  TCanvas* c1 = new TCanvas("c1","c1");
  effVwVlenx->Divide(hitsVwVlenx,cellsVwVlenx);
  effVwVlenx->Draw("colz");
  effVwVlenx->GetZaxis()->SetRangeUser(0.5,1.0);

  TCanvas* c2 = new TCanvas("c2","c2");
  effVwVleny->Divide(hitsVwVleny,cellsVwVleny);
  effVwVleny->Draw("colz");
  effVwVleny->GetZaxis()->SetRangeUser(0.5,1.0);

  TCanvas* c3 = new TCanvas("c3","c3");
  effVwVlenratio->Divide(effVwVleny,effVwVlenx);
  effVwVlenratio->Draw("colz");
  effVwVlenratio->GetZaxis()->SetRangeUser(0.8,1.2);


  double nxcell = 0.0;
  double nycell = 0.0;
  
  double nxcellcut = 0.0;
  double nycellcut = 0.0;
  
  // loop over every cell in the first diblock
  for(int icell = 1; icell <= ncellbins; ++icell){
    // convert this to a plane and cell number
    int cell = (icell-1)%384;
    int plane = (icell-1-cell)/384;
    int modcell = cell%32;
//     int module = (cell-modcell)/32;
    
    int view = plane%2; // view = 0 y, view = 1 x
    // need to switch the module in cell for x view to keep consistent cell numbering within module between views
    if(view){ modcell = 31-modcell; }

    int celltot = cellsVwcell->ProjectionY("",icell,icell,"")->GetEntries();
    // If never get any tri-cells on this cell skip it    
    if(celltot == 0){ continue; }
    if(plane > (int)plmax){ continue; }

    if(view){ ++nxcell; }
    else{ ++nycell; }
    if(celltot > nmincell){
      if(view){ ++nxcellcut; }
      else{ ++nycellcut; }
    }
    // loop over every w bin
    for(int iw = 1; iw <= nwbins; ++iw){
      // get the efficiency at this cell/w, eff = k/n
      double k = hitsVwcell->GetBinContent(icell,iw);
      double n = cellsVwcell->GetBinContent(icell,iw);
      double ef = k/n;
      if(n == 0){ ef = 0; }
      // calculate the variance of this
      double var = (k+1)*(k+2)/((n+2)*(n+3)) - (k+1)*(k+1)/((n+2)*(n+2));

      double w = hitsVwcell->GetYaxis()->GetBinCenter(iw);
      
      eff->SetBinContent(iw,ef);
      if(view){
        xeffVw->SetBinContent(iw,(ef+xeffVw->GetBinContent(iw)));
        effVwx->Fill(w,ef);
        if(celltot > nmincell){
          xeffVwcut->SetBinContent(iw,(ef+xeffVwcut->GetBinContent(iw)));
          effVwxcut->Fill(w,ef);
        }
      }
      else{
        yeffVw->SetBinContent(iw,(ef+yeffVw->GetBinContent(iw)));
        effVwy->Fill(w,ef);
        if(celltot > nmincell){
          yeffVwcut->SetBinContent(iw,(ef+yeffVwcut->GetBinContent(iw)));
          effVwycut->Fill(w,ef);
        }
      }
      
      eff->SetBinError(iw,sqrt(var));      
    } // end loop over w bins
    
    // now get the average efficiency of this cell over different w ranges
    double lowedgeeff = eff->Integral(1,lowwbin)/((double)(lowwbin - 1 + 1));
    double mideff = eff->Integral(lowwbin+1,highwbin-1)/((double)(highwbin-1 - (lowwbin+1) + 1));
    double highedgeeff = eff->Integral(highwbin,nwbins)/((double)(nwbins - (highwbin) + 1));

    int planebin = effslowCellVPlanex->GetXaxis()->FindFixBin(plane);
    int cellbin = effslowCellVPlanex->GetYaxis()->FindFixBin(cell);

    if(view){
      effslowCellVPlanex->SetBinContent(planebin,cellbin,lowedgeeff);
      effslowx->Fill(lowedgeeff);
      effsmidCellVPlanex->SetBinContent(planebin,cellbin,mideff);
      effsmidx->Fill(mideff);
      effshighCellVPlanex->SetBinContent(planebin,cellbin,highedgeeff);
      effshighx->Fill(highedgeeff);
    }
    else{
      effslowCellVPlaney->SetBinContent(planebin,cellbin,lowedgeeff);
      effslowy->Fill(lowedgeeff);
      effsmidCellVPlaney->SetBinContent(planebin,cellbin,mideff);
      effsmidy->Fill(mideff);
      effshighCellVPlaney->SetBinContent(planebin,cellbin,highedgeeff);
      effshighy->Fill(highedgeeff);
    }
        
  }
  
  
  // scale all the histograms based on the number of entries
  xeffVw->Scale(1.0/nxcell);
  yeffVw->Scale(1.0/nycell);
  xeffVwcut->Scale(1.0/nxcellcut);
  yeffVwcut->Scale(1.0/nycellcut);
  
    
  
  // set errors  
  for(int iw = 1; iw <= xeffVw->GetNbinsX(); ++iw){
    double err = xeffVw->GetBinError(iw);
    xeffVw->SetBinError(iw,(sqrt(err)/nxcell));
  }
  for(int iw = 1; iw <= yeffVw->GetNbinsX(); ++iw){
    double err = yeffVw->GetBinError(iw);
    yeffVw->SetBinError(iw,(sqrt(err)/nycell));
  }
  TCanvas* c4 = new TCanvas("c4","c4");
  xeffVw->Draw();
  yeffVw->Draw("same");
  xeffVw->GetYaxis()->SetRangeUser(0.0,1.1);
  xeffVw->SetLineColor(kBlue);
  yeffVw->SetLineColor(kRed);
  xeffVw->SetMarkerColor(kBlue);
  yeffVw->SetMarkerColor(kRed);
  TLegend* leg4 = new TLegend(0.55,0.15,0.85,0.35);
  leg4->AddEntry(xeffVw,"Vertical Cells","LE");
  leg4->AddEntry(yeffVw,"Horizontal Cells","LE");
  leg4->Draw("same");

  TCanvas* c5 = new TCanvas("c5","c5");
  TProfile* effVwprofx = effVwx->ProfileX("effVwprofx",1,-1,"s");
  effVwprofx->GetYaxis()->SetRangeUser(0.0,1.1);
  effVwprofx->GetYaxis()->SetTitle("Efficiency");
  effVwprofx->Draw();
  TProfile* effVwprofy = effVwy->ProfileX("effVwprofy",1,-1,"s");
  effVwprofy->Draw("same");
  effVwprofx->SetLineColor(kBlue);
  effVwprofx->SetMarkerColor(kBlue);
  effVwprofy->SetLineColor(kRed);
  effVwprofy->SetMarkerColor(kRed);
  TLegend* leg5 = new TLegend(0.55,0.15,0.85,0.35);
  leg5->AddEntry(effVwprofx,"Vertical Cells","LE");
  leg5->AddEntry(effVwprofy,"Horizontal Cells","LE");
  leg5->Draw("same");
    
  
  for(int iw = 1; iw <= xeffVwcut->GetNbinsX(); ++iw){
    double err = xeffVwcut->GetBinError(iw);
    xeffVwcut->SetBinError(iw,(sqrt(err)/nxcellcut));
  }
  for(int iw = 1; iw <= yeffVwcut->GetNbinsX(); ++iw){
    double err = yeffVwcut->GetBinError(iw);
    yeffVwcut->SetBinError(iw,(sqrt(err)/nycellcut));
  }
  TCanvas* c6 = new TCanvas("c6","c6");
  xeffVwcut->Draw();
  yeffVwcut->Draw("same");
  xeffVwcut->GetYaxis()->SetRangeUser(0.0,1.1);
  xeffVwcut->SetLineColor(kBlue);
  yeffVwcut->SetLineColor(kRed);
  xeffVwcut->SetMarkerColor(kBlue);
  yeffVwcut->SetMarkerColor(kRed);
  TLegend* leg6 = new TLegend(0.55,0.15,0.85,0.35);
  leg6->AddEntry(xeffVwcut,"Vertical Cells","LE");
  leg6->AddEntry(yeffVwcut,"Horizontal Cells","LE");
  leg6->Draw("same");
  
  TCanvas* c7 = new TCanvas("c7","c7");
  TProfile* effVwprofxcut = effVwxcut->ProfileX("effVwprofxcut",1,-1,"s");
  effVwprofxcut->GetYaxis()->SetRangeUser(0.0,1.1);
  effVwprofxcut->GetYaxis()->SetTitle("Efficiency");
  effVwprofxcut->Draw();
  TProfile* effVwprofycut = effVwycut->ProfileX("effVwprofycut",1,-1,"s");
  effVwprofycut->Draw("same");
  effVwprofxcut->SetLineColor(kBlue);
  effVwprofxcut->SetMarkerColor(kBlue);
  effVwprofycut->SetLineColor(kRed);
  effVwprofycut->SetMarkerColor(kRed);
  TLegend* leg7 = new TLegend(0.55,0.15,0.85,0.35);
  leg7->AddEntry(effVwprofxcut,"Vertical Cells","LE");
  leg7->AddEntry(effVwprofycut,"Horizontal Cells","LE");
  leg7->Draw("same");
        
  TCanvas* c8 = new TCanvas("c8","c8");
  effslowCellVPlanex->Draw("colz");
  effslowCellVPlanex->GetZaxis()->SetRangeUser(0.0,1.0);
  
  TCanvas* c9 = new TCanvas("c9","c9");
  effslowCellVPlaney->Draw("colz");
  effslowCellVPlaney->GetZaxis()->SetRangeUser(0.0,1.0);
  
  TCanvas* c10 = new TCanvas("c10","c10");
  effslowx->Draw();
  effslowy->Draw("same");
  effslowx->SetLineColor(kBlue);
  effslowy->SetLineColor(kRed);
  TLegend* leg10 = new TLegend(0.15,0.65,0.35,0.85);
  leg10->AddEntry(effslowx,"X Cells","LE");
  leg10->AddEntry(effslowy,"Y Cells","LE");
  leg10->Draw("same");
  
  TCanvas* c11 = new TCanvas("c11","c11");
  effsmidCellVPlanex->Draw("colz");
  effsmidCellVPlanex->GetZaxis()->SetRangeUser(0.0,1.0);
  
  TCanvas* c12 = new TCanvas("c12","c12");
  effsmidCellVPlaney->Draw("colz");
  effsmidCellVPlaney->GetZaxis()->SetRangeUser(0.0,1.0);
  
  TCanvas* c13 = new TCanvas("c13","c13");
  effsmidx->Draw();
  effsmidy->Draw("same");
  effsmidx->SetLineColor(kBlue);
  effsmidy->SetLineColor(kRed);
  TLegend* leg13 = new TLegend(0.15,0.65,0.35,0.85);
  leg13->AddEntry(effsmidx,"X Cells","LE");
  leg13->AddEntry(effsmidy,"Y Cells","LE");
  leg13->Draw("same");
  
  TCanvas* c14 = new TCanvas("c14","c14");
  effshighCellVPlanex->Draw("colz");
  effshighCellVPlanex->GetZaxis()->SetRangeUser(0.0,1.0);
  
  TCanvas* c15 = new TCanvas("c15","c15");
  effshighCellVPlaney->Draw("colz");
  effshighCellVPlaney->GetZaxis()->SetRangeUser(0.0,1.0);
  
  TCanvas* c16 = new TCanvas("c16","c16");
  effshighx->Draw();
  effshighy->Draw("same");
  effshighx->SetLineColor(kBlue);
  effshighy->SetLineColor(kRed);
  TLegend* leg16 = new TLegend(0.15,0.65,0.35,0.85);
  leg16->AddEntry(effshighx,"X Cells","LE");
  leg16->AddEntry(effshighy,"Y Cells","LE");
  leg16->Draw("same");
      
}