MakeMEUPerFEB.C

Go to the documentation of this file.

#include "CalibAnaPlot.h"
#include "CalibAnaTypes.h"
#include "TexBuilder.h"
#include <iomanip>
#include <numeric>

class MEUPerFEB : public CalibAnaPlot{

  public:

    MEUPerFEB(std::string outpath);
    // vvv CalibAnaPlot implimentations
    void SchedulePlots();
    void ScheduleAllDataPlots();
    void ScheduleSamples();
    void ScheduleHitVars();
    void FillHitPlots();
    // ^^^ CalibAnaPlot implimentations
    void MakeTex();
    void InitBadFeb();
    bool BadFeb(int plane, int cell);
    void DriftTex(std::string plotdir);

  

  private:
    std::vector< std::string > fEpochList;
    TexBuilder * fTex;
    std::map< int, std::map < int, bool > > fBadFeb;//[plane][module]
};


MEUPerFEB::MEUPerFEB(std::string outpath)
: CalibAnaPlot(outpath)
{
}

void MakeMEUPerFEB(std::string outpath,
                   std::string sample,
                   bool justdraw=false,
                   unsigned int stride=0,
                   unsigned int limit=0)
{
  MEUPerFEB meuperapd(outpath);
  meuperapd.Initialize();
  if(stride) meuperapd.SetTreeLoopStride(stride);
  if(limit) meuperapd.SetTreeLoopLimit(limit);

  if(justdraw){
    meuperapd.MakeTex();
  } else {
    if(sample=="all"){
      for( auto& s : meuperapd.GetSamples() ){
        meuperapd.ProcessTrees(s.first);
      } 
      meuperapd.MakeTex(); 
    } else {
      meuperapd.ProcessTrees(sample);
    }
  }
}

//^~^~^~~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void MEUPerFEB::MakeTex()
{
  InitBadFeb();

  fTex = new TexBuilder( fOutPath+"/MEUPerCellGroupAndFEB.tex", 
                         "Drift by FEB", 
                         "Tyler Alion, Alexander Booth, Yibing Zhang",
                         "Miniproduction 5 Data vs. MC. Tricell hits from Stopping Cosmics populate all plots, with Absolute Calibration cuts applied. Showing PE/cm and PECorr/cm broken down by cell group and view then further by FEB.");
  
  std::string plotdir = "nd_all_samples";

  // There are 12 denominator plots
  std::vector< std::string > views = {"x", "y"};
  std::vector< std::string > mods  = {"_0-31", "_32-63", "_64-95"};
  std::vector< std::string > yaxes = {"pecm", "pecorrcm"};
  for(auto& view : views){
    size_t firstplane = 1;
    if(view=="y") firstplane = 0;
    for(auto& modstr : mods){
      for(auto& yax : yaxes){

        std::vector< std::string > drawplots;
        size_t np = fDet->NPlanes();
        for(size_t p=firstplane; p<np; p+=2)
          drawplots.push_back(yax+"_time_"+view+"_"+std::to_string(p)+modstr);
        drawplots.push_back(yax+"_time_"+view+modstr);
        drawplots.push_back(yax+"_time_"+view+modstr+"_febcut");
        
        Plot denom_plot(fOutPath+"/nd_all_data.root",         // data over all time
                        yax+"_time_"+view+modstr+"_febcut",   // mean for whole view-module group
                        "Mean: "+view+modstr);

        for(auto& plotname : drawplots){        
          std::vector< Plot > plotstack;
          for(auto& samp : fSamples)
            plotstack.push_back( Plot( fOutPath+"/"+samp.first+".root", 
                                       plotname,
                                       LegendLabel(samp.first) ) );
          
          Draw( plotdir, plotstack, denom_plot, "" );
        }
      }
    }
  }

  

  DriftTex(plotdir);

  delete fTex;
}

//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void MEUPerFEB::DriftTex(std::string plotdir)
{
  fTex->SectionNoNum("By View and Cell Group");
  fTex->AddFigureRow({plotdir+"/legend"}, {""}, {""});
  fTex->AddFigureRow({plotdir+"/pecm_time_x_0-31",  plotdir+"/pecm_time_y_0-31" }, {"", ""}, {"", ""}); 
  fTex->AddFigureRow({plotdir+"/pecm_time_x_32-63", plotdir+"/pecm_time_y_32-63"}, {"", ""}, {"", ""}); 
  fTex->AddFigureRow({plotdir+"/pecm_time_x_64-95", plotdir+"/pecm_time_y_64-95"}, {"", ""}, {"", ""}); 
  fTex->NewPage();
  fTex->AddFigureRow({plotdir+"/legend"}, {""}, {""});
  fTex->AddFigureRow({plotdir+"/pecorrcm_time_x_0-31",  plotdir+"/pecorrcm_time_y_0-31" }, {"", ""}, {"", ""}); 
  fTex->AddFigureRow({plotdir+"/pecorrcm_time_x_32-63", plotdir+"/pecorrcm_time_y_32-63"}, {"", ""}, {"", ""}); 
  fTex->AddFigureRow({plotdir+"/pecorrcm_time_x_64-95", plotdir+"/pecorrcm_time_y_64-95"}, {"", ""}, {"", ""}); 
  fTex->NewPage();

  fTex->SectionNoNum("MEU with FEB Cut");
  fTex->AddFigureRow({plotdir+"/legend"}, {""}, {""});
  fTex->AddFigureRow({plotdir+"/pecm_time_x_0-31_febcut",  plotdir+"/pecm_time_y_0-31_febcut" }, {"", ""}, {"", ""}); 
  fTex->AddFigureRow({plotdir+"/pecm_time_x_32-63_febcut", plotdir+"/pecm_time_y_32-63_febcut"}, {"", ""}, {"", ""}); 
  fTex->AddFigureRow({plotdir+"/pecm_time_x_64-95_febcut", plotdir+"/pecm_time_y_64-95_febcut"}, {"", ""}, {"", ""}); 
  fTex->NewPage();
  fTex->AddFigureRow({plotdir+"/legend"}, {""}, {""});
  fTex->AddFigureRow({plotdir+"/pecorrcm_time_x_0-31_febcut",  plotdir+"/pecorrcm_time_y_0-31_febcut" }, {"", ""}, {"", ""}); 
  fTex->AddFigureRow({plotdir+"/pecorrcm_time_x_32-63_febcut", plotdir+"/pecorrcm_time_y_32-63_febcut"}, {"", ""}, {"", ""}); 
  fTex->AddFigureRow({plotdir+"/pecorrcm_time_x_64-95_febcut", plotdir+"/pecorrcm_time_y_64-95_febcut"}, {"", ""}, {"", ""}); 
  fTex->NewPage();


  unsigned int nPlanes = fDet->NPlanes();
  for(unsigned int p = 0; p < nPlanes; p++)
  {
    fTex->SectionNoNum(Form("Plane %i", p));
    std::string pstr = (p % 2 == 0 ? "y_" : "x_") + std::to_string(p);
    fTex->AddFigureRow({plotdir+"/legend"}, {""}, {""});
    fTex->AddFigureRow({plotdir+"/pecm_time_"+pstr+"_0-31",  plotdir+"/pecorrcm_time_"+pstr+"_0-31" }, {"", ""}, {"", ""}); 
    fTex->AddFigureRow({plotdir+"/pecm_time_"+pstr+"_32-63", plotdir+"/pecorrcm_time_"+pstr+"_32-63"}, {"", ""}, {"", ""}); 
    fTex->AddFigureRow({plotdir+"/pecm_time_"+pstr+"_64-95", plotdir+"/pecorrcm_time_"+pstr+"_64-95"}, {"", ""}, {"", ""}); 
    if(p!=nPlanes-1) fTex->NewPage();
  }
}

//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void MEUPerFEB::FillHitPlots()
{
  // The filled branches are protected members of
  // CalibAnaPlot, so can be used here

  bool abscal_hit = ( w > -100 && w < 100 && 
                      path > 0. &&
                      pecorr > 0. && pecorr/path < 100. &&
                      pe > 0. &&
                      cmFromEnd > 100 && cmFromEnd < 200 );
  if( !abscal_hit ) return;
  if(  plane > 191) return;

  std::string vstr, pstr;
  if(view==0){ vstr="x"; }
  if(view==1){ vstr="y"; }  
  pstr = std::to_string(plane);

  std::string end      = vstr + "_" + pstr + "_";
  std::string endGroup = vstr + "_";

  float pecm     = pe/path;
  float pecorrcm = pecorr/path;
  if(cell <= 31){
    end      += "0-31";
    endGroup += "0-31";
  }
  else if(cell > 31 && cell <= 63){
    end      += "32-63";
    endGroup += "32-63";
  }
  else{
    end      += "64-95";
    endGroup += "64-95";
  }
  fPlots["pecm_time_"+end]    ->Fill(evt_time,pecm);
  fPlots["pecorrcm_time_"+end]->Fill(evt_time,pecorrcm);
  fPlots["pecm_time_"+endGroup]    ->Fill(evt_time,pecm);
  fPlots["pecorrcm_time_"+endGroup]->Fill(evt_time,pecorrcm);

  fAllDataPlots["pecm_time_"+endGroup]    ->Fill(evt_time,pecm);
  fAllDataPlots["pecorrcm_time_"+endGroup]->Fill(evt_time,pecorrcm);
  
  if( !BadFeb(plane,cell) ){
    fPlots[        "pecm_time_"+endGroup+"_febcut"     ]->Fill(evt_time,pecm);
    fPlots[        "pecorrcm_time_"+endGroup+"_febcut" ]->Fill(evt_time,pecorrcm);      
    fAllDataPlots[ "pecm_time_"+endGroup+"_febcut"     ]->Fill(evt_time,pecm);
    fAllDataPlots[ "pecorrcm_time_"+endGroup+"_febcut" ]->Fill(evt_time,pecorrcm);
  }
}

//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
bool MEUPerFEB::BadFeb(int plane, int cell)
{
  // If no plane entry, it is _not_ a bad FEB
  if( !fBadFeb.count(plane) ) return false;
  // If there is a plane entry...
  return fBadFeb[plane].count( fDet->CellToModule(cell) );
}

//^~^~^~~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void MEUPerFEB::InitBadFeb()
{
  // From a hand scan.

  // Eventually use algorith which compares Mean&RMS of time bin
  //  in a given FEB to it's view-module average 

  // Y view
  std::vector< int > planes = { 2, 10, 14, 16, 24, 38, 46, 50, 66, 82, 88, 90, 94, 114, 140, 150, 152, 184 };
  for(auto& p : planes)
    fBadFeb[p][1]=true;
  planes = { 12, 20, 46, 110, 124, 130, 166 };
  for(auto& p : planes)
    fBadFeb[p][2]=true;

  // X view
  planes = { 5, 63, 75, 77, 93, 95, 135, 149 };
  for(auto& p : planes)
    fBadFeb[p][1]=true;
  planes = { 17, 31, 35, 69, 91, 95, 129, 175 };
  for(auto& p : planes)
    fBadFeb[p][2]=true;

}

//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void MEUPerFEB::ScheduleHitVars()
{
  fHitVars="w, view, plane, cell, pe, pecorr, path, evt_time, cmFromEnd,";
}

//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void MEUPerFEB::SchedulePlots()
{
  int nPlanes = fDet->NPlanes();

  std::vector< PlotAxis > xaxes = { 
    {"time",  "Unix Time (s)",  100, 14e8, 16e8}
  };

  std::vector< ContentAxis > yaxes = { 
    {"pecm",     "PE/cm"},
    {"pecorrcm", "PECorr/cm"},
  };

  std::vector< std::string > cells = {"0-31", "32-63", "64-95"};
  std::vector<int> planes(nPlanes);
  std::iota(planes.begin(), planes.end(), 0);
  for(auto& p : planes){
    std::string vstr = (p % 2 == 0 ? "Y" : "X");
    std::string v    = (p % 2 == 0 ? "y" : "x");
    for(auto& c : cells){
      for(auto& xax : xaxes){
        for(auto& yax : yaxes){

          std::string name = yax.label +"_"+ xax.label +"_"+ v + "_" + std::to_string(p) + "_" + c;
          //std::cout << "Schedule: " << name << std::endl;

          fPlots[name] = new TProfile( name.c_str(),
              (vstr+" View, Plane "+std::to_string(p)+", Cells "+c+";" + xax.title +";"+ yax.title).c_str(),
              xax.nbins, xax.min, xax.max );

          if(p < 2){
            name = yax.label +"_"+ xax.label +"_"+ v + "_" + c;
            fPlots[name] = new TProfile( name.c_str(),
                (vstr+" View, Cells "+c+";" + xax.title +";"+ yax.title).c_str(),
                xax.nbins, xax.min, xax.max );
            name = name + "_febcut";
            fPlots[name] = new TProfile( name.c_str(),
                (vstr+" View, Cells "+c+" + FEB Cut;" + xax.title +";"+ yax.title).c_str(),
                xax.nbins, xax.min, xax.max );
          } // by module
        } // yaxes
      } // xaxes
    } // cell groups
  } // planes
} // SchedulePlots

//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void MEUPerFEB::ScheduleAllDataPlots()
{

  std::vector< PlotAxis > xaxes = { 
    {"time",  "Unix Time (s)",  100, 14e8, 16e8}
  };

  std::vector< ContentAxis > yaxes = { 
    {"pecm",     "PE/cm"},
    {"pecorrcm", "PECorr/cm"},
  };

  std::vector< std::string > cells = {"0-31", "32-63", "64-95"};
  std::vector< std::string > views = {"x", "y"};
  for(auto& v : views){
    for(auto& c : cells){
      for(auto& xax : xaxes){
        for(auto& yax : yaxes){
          std::string name = yax.label +"_"+ xax.label +"_"+v+"_" + c;
          std::string vCap = (v=="x") ? "X" : "Y";
          fAllDataPlots[name] = new TProfile( name.c_str(),
                                              (vCap+" View, Cells "+c+";" + xax.title +";"+ yax.title).c_str(),
                                              xax.nbins, xax.min, xax.max );
          name = name + "_febcut";
          fAllDataPlots[name] = new TProfile( name.c_str(),
                                              (vCap+" View, Cells "+c+" + FEB Cut;" + xax.title +";"+ yax.title).c_str(),
                                              xax.nbins, xax.min, xax.max );

        }
      }
    }
  }
}

//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void MEUPerFEB::ScheduleSamples()
{
  std::vector< std::string > epoch_list = {
    //"4c"};    
    "2a", "2b","2c","2d","2e",
    "3b","3c","3d","3e",
    "4a", "4b.1", "4b.2", "4c", // split 4b up manually
    "5a","5b",
    "6a","6b",
    "7a","7b","7c","7d",
    "8b" };

  fEpochList = epoch_list;

  for( auto& ep : epoch_list ){
    if(ep=="4b.1" || ep=="4b.2") 
      fSamples["nd_data_ep"+ep]="/nova/ana/users/talion/miniprod5/pcliststop_miniprod5.b/nd_data_epoch4b/calibana/merged.root";
    else 
      fSamples["nd_data_ep"+ep]="/nova/ana/users/talion/miniprod5/pcliststop_miniprod5.b/nd_data_epoch"+ep+"/calibana/merged.root";
  }
  
  //fSamples["nd_mc"]="/nova/ana/users/abooth/miniprod5.b/nd_mc/merged.root";
}