NCExtractionCuts.cxx

Go to the documentation of this file.

//
//
//NCExtractionCuts.cxx
//
//Module for analyzing beam data for nc analysis
//
//T. Osiecki 9/2006

#include "NCUtils/Extraction/NCExtractionCuts.h"

#include "MessageService/MsgService.h"
#include "Conventions/Detector.h"
#include "Conventions/SimFlag.h"
#include "NCUtils/Cuts/NCAnalysisCutsNC.h"

#include "TDirectory.h"

#include <cassert>

#include "NCUtils/Extraction/MicroDSTMaker.h"
REGISTER_NCEXTRACTION(NCExtractionCuts, Cuts)

static const int kNumPDFBaseNamesDP = 3;

ClassImp(NCExtractionCuts)

CVSID("$Id: NCExtractionCuts.cxx,v 1.10 2009/11/23 16:55:22 rodriges Exp $");

//......................................................................
NCExtractionCuts::NCExtractionCuts(NCAnalysisCuts* cuts,
                                   const Registry& r) :
  NCExtraction(cuts, r),
  fPDFsFilled(false),
  fBeamType(BeamType::kUnknown)

{
  MSG("JobC", Msg::kInfo) << "NCExtractionCuts::Constructor 2" << endl;


  const char* tmps;
  if (r.Get("BeamType", tmps)) fBeamType=BeamType::TagToEnum(tmps);
}

//---------------------------------------------------------------------------
void NCExtractionCuts::ReadPDFs(NCEventInfo& evtInfo)
{
  MSG("NCExtractionCuts", Msg::kInfo) << "ReadPDFs" << endl;

  TDirectory *savedir = gDirectory;

  // Try $SRT_LOCAL, $SRT_PRIVATE_CONTEXT then $SRT_PUBLIC_CONTEXT
  TString basedir(getenv("SRT_LOCAL"));
  if(basedir=="") basedir=getenv("SRT_PRIVATE_CONTEXT");
  // SRT_PRIVATE_CONTEXT gets set to '.' instead of empty, for some reason
  if(basedir==".") basedir=getenv("SRT_PUBLIC_CONTEXT");

  TString pdfPath = basedir+"/NCUtils/data/";
  
  Detector::Detector_t detType = (Detector::Detector_t)evtInfo.header->detector;
  
  TString pdffile = WhichCCPIDFiles(detType, fReleaseType, fBeamType);
  assert(pdffile!="");
  
  pdfPath += pdffile;
  MSG("NCExtractionCuts", Msg::kInfo)
    << "Reading PDFS from file " << pdfPath << endl;
  savedir->cd();
  
  TFile pdfFile(pdfPath);
  
  //clear the existing vectors
  fPDFs.clear();
  const char* pdfNamesInFile[kNumPDFBaseNamesDP] = {"Event length - %s",
                                                    "Track ph frac - %s",
                                                    "ph per plane (%s)"};
  vector<TString> ncccNames;
  ncccNames.resize(2);
  ncccNames[NCType::kNC]="nc";
  ncccNames[NCType::kCC]="cc";

  for(int i=0; i<kNumPDFBaseNamesDP; ++i){
    for(int evtType=0; evtType<2; ++evtType){
      TString histname=Form(pdfNamesInFile[i], ncccNames[evtType].Data());
      TH1F* hist=(TH1F*)pdfFile.Get(histname);
      if(!hist){ 
        MSG("NCExtractionCuts", Msg::kFatal) 
          << "Can't find histogram with name " << histname 
          << " in file." << endl;
        // Dummy message to make sure it aborts
        MSG("NCExtractionCuts", Msg::kFatal) << "Goodbye" << endl;
      }
      hist->SetDirectory(0);
      fPDFs[(NCType::EEventType)evtType].push_back(hist);
    }
  }

  savedir->cd();

  fPDFsFilled = true;
}

//----------------------------------------------------------------------
float NCExtractionCuts::FindPDFProbability(TH1F *pdf, float pdfValue)
{
  int bin = 0;
  float binContent = 0.;

  bin = pdf->FindBin(pdfValue);
  binContent = pdf->GetBinContent(bin);

  if(binContent>0. && pdf->Integral()>0.) return binContent/pdf->Integral();
  else                                    return 0;
}

//-----------------------------------------------------------------------
double NCExtractionCuts::GetIdProbability(NCEventInfo& evtInfo, int)
{
  MAXMSG("NCExtractionCuts", Msg::kDebug,20) << "GetIdProbabilty" << endl;

  bool passesNCCuts=PassesNCCuts(evtInfo);
  bool passesCCCuts=PassesCCCuts(evtInfo);

  // Mini truth table for return value:
  //     |  CC  | !CC
  //  ----------------
  //  NC |  -1  |   0 
  // !NC |   1  |  -2

  if     ( passesCCCuts && !passesNCCuts) return 1;
  else if(!passesCCCuts &&  passesNCCuts) return 0;
  else if( passesCCCuts &&  passesNCCuts) return -1;
  else if(!passesCCCuts && !passesNCCuts) return -2;

  assert(0 && "This should never happen");
}

//-----------------------------------------------------------------------
bool NCExtractionCuts::PassesNCCuts(NCEventInfo& evtInfo)
{
  const Detector::Detector_t det = (Detector::Detector_t)evtInfo.header->detector;

  int tracks = evtInfo.event->tracks;

  int planes = det==Detector::kNear ? evtInfo.event->lengthInPlanes :
                                      evtInfo.event->planes;

  int trackExtension = evtInfo.shower->planes - evtInfo.track->planes;

  //NC Algorithm:

  // 1. Long events are definitely not NC
  if(planes >= 47) return false;
  else{
    // 2. Events without tracks definitely *are* NC
    if(tracks == 0) {
      return true;
    } else {
      // 3. Short events with tracks need to pass the track extension cut
      if(trackExtension > -6) {
        return true;
      } else {
        return false;
      }
    }
  }

}

//-----------------------------------------------------------------------
bool NCExtractionCuts::PassesCCCuts(NCEventInfo& evtInfo)
{
  const Detector::Detector_t det = (Detector::Detector_t)evtInfo.header->detector;

  // David uses end-beg for the event length in the near,
  //not the number of hit planes
  int dpPlanes = det==Detector::kNear ? evtInfo.event->lengthInPlanes :
                                        evtInfo.event->eventSummaryPlanes;
  
  float trackSignalUseFraction = evtInfo.track->pulseHeight;
  if(evtInfo.event->pulseHeight > 0.)
    trackSignalUseFraction /= evtInfo.event->pulseHeight;

  float trackSignalPerPlane = evtInfo.track->pulseHeight;
  if(evtInfo.track->planes > 0) trackSignalPerPlane /= evtInfo.track->planes;

  vector<float> dpInputs;

  dpInputs.push_back(dpPlanes);
  dpInputs.push_back(trackSignalUseFraction);
  dpInputs.push_back(trackSignalPerPlane);

  float ncLikelihood = 1.;
  float ccLikelihood = 1.;

  for(int i = 0; i < (int)dpInputs.size(); ++i){
    //check to see if we should use the current distribution
    float ncProbability = FindPDFProbability(fPDFs[NCType::kNC][i],
                                             dpInputs[i]);
    float ccProbability = FindPDFProbability(fPDFs[NCType::kCC][i],
                                             dpInputs[i]);

    if(ncProbability == 0.) ncProbability = 1e-4;
    if(ccProbability == 0.) ccProbability = 1e-4;

    ncLikelihood *= ncProbability;
    ccLikelihood *= ccProbability;
  }//end over pdfs



  double probPID = TMath::Sqrt(-TMath::Log(ccLikelihood)) - 
                   TMath::Sqrt(-TMath::Log(ncLikelihood));
  probPID *= -1.;

  const float cc_cut = det == Detector::kNear ? -0.1 : -0.2;

  int tracks = evtInfo.event->tracks;

  bool passesCCCuts = probPID > cc_cut && tracks > 0;
  // In the near detector we additionally require to pass the track
  // fit pass/reclamation cut
  if(det == Detector::kNear) 
    passesCCCuts = passesCCCuts && PassesTrackReclamation(evtInfo);

  return passesCCCuts;
}

//-----------------------------------------------------------------------
bool NCExtractionCuts::PassesTrackReclamation(NCEventInfo& evtInfo)
{
  return evtInfo.track->passedFit>0 ||
    (TMath::Abs(evtInfo.track->begPlaneU - evtInfo.track->begPlaneV) < 6
     && TMath::Abs(evtInfo.track->endPlaneU - evtInfo.track->endPlaneV) < 41
     && evtInfo.track->endPlane < 270);
}

//-----------------------------------------------------------------------
void NCExtractionCuts::FillAnalysisInfo(NCEventInfo& evtInfo, int beamType)
{

  if(!fPDFsFilled) ReadPDFs(evtInfo);

  ANtpAnalysisInfo* ana=evtInfo.analysis;
  //double separationParameter=GetIdProbability(evtInfo, beamType);
  double separationCut=GetCutPosition();

  ana->Reset();

  ana->separationParameterCut = separationCut;
  ana->separationParameter = GetIdProbability(evtInfo, beamType);

  // Events that pass the NC algorithm are NC...
  ana->isNC = PassesNCCuts(evtInfo);
  // ...of those that aren't, only the ones passing the CC cuts are CC
  ana->isCC = (!ana->isNC) && PassesCCCuts(evtInfo);

  if(ana->isCC) {
    ana->pass = (int)( evtInfo.reco->trackMomentum <= 0 && 
                       fCuts->IsGoodBeamEvent() && 
                       fCuts->InBeamFiducialVolume() );
  } else if(ana->isNC) {
    ana->pass = (int)fCuts->InBeamFiducialVolume();
  }
  else ana->pass=0; // neither NC nor CC
}

//-----------------------------------------------------------------------
TString NCExtractionCuts::WhichCCPIDFiles(Detector::Detector_t detType,
                                          ReleaseType::Release_t mcType,
                                          BeamType::BeamType_t beamType)
{

  if(ReleaseType::IsDogwood(mcType)){
    MAXMSG("NCExtractionCuts", Msg::kError,5)                   << endl
      << "****************************************************" << endl
      << "* USING CEDAR CC PID FILES ON DOGWOOD. FIX THIS    *" << endl
      << "* BEFORE DOING A REAL ANALYSIS                     *" << endl
      << "****************************************************" << endl;
  }

  assert(mcType != ReleaseType::kUnknown);

  assert(!ReleaseType::IsCarrot(mcType));

  if(detType == Detector::kFar) {
     if(ReleaseType::IsDaikon(mcType)) {
       switch (beamType){
       case BeamType::kL010z185i:
       case BeamType::kL010z185i_i124:
       case BeamType::kL010z185i_i191:
       case BeamType::kL010z185i_i213:
       case BeamType::kL010z185i_i224:
       case BeamType::kL010z185i_i232:
       case BeamType::kL010z185i_i243:
       case BeamType::kL010z185i_i257:
       case BeamType::kL010z185i_i282:
       case BeamType::kL010z185i_i303:
       case BeamType::kL010z185i_i324: 
       case BeamType::kUnknown:
         return "dp_pdf_far_le_cedar_daikon.root";
       default:
         return "";
       }
     }
     return "";
  }

  if(detType == Detector::kNear){
    if(ReleaseType::IsDaikon(mcType)){
      switch (beamType){
      case BeamType::kL000z200i:
      case BeamType::kL010z185i:
      case BeamType::kL010z185i_i124:
      case BeamType::kL010z185i_i191:
      case BeamType::kL010z185i_i213:
      case BeamType::kL010z185i_i224:
      case BeamType::kL010z185i_i232:
      case BeamType::kL010z185i_i243:
      case BeamType::kL010z185i_i257:
      case BeamType::kL010z185i_i282:
      case BeamType::kL010z185i_i303:
      case BeamType::kL010z185i_i324:
      case BeamType::kL050z200i:
      case BeamType::kL200z200i:
      case BeamType::kL010z000i:
      case BeamType::kL010z170i:
      case BeamType::kL010z200i:
      case BeamType::kL010z185i_lowintensity:
      case BeamType::kUnknown:
        return "dp_pdf_near_L010z185i_cedar_daikon.root";

      case BeamType::kL100z200i:
        return "dp_pdf_near_L100z200i_cedar_daikon.root";

      case BeamType::kL250z200i:
        return "dp_pdf_near_L250z200i_cedar_daikon.root";

      case BeamType::kL150z200i:
        return "dp_pdf_near_L150z200i_cedar_daikon.root";

      default:
        return "";
      }
    }
    return "";
  }

  return "";
}

---