ValVtxModule.cxx

Go to the documentation of this file.

// $Id: ValVtxModule.cxx,v 1.4 2008/07/11 19:19:48 boehm Exp $
//
// FILL_IN: [Document your code!!]
//
// vahle@hep.ucl.ac.uk
#include <dirent.h>
#include "TChain.h"
#include "TFile.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "VertexFinder/Module/VtxRecord.h"
#include "VertexFinder/Module/ValVtxModule.h"
#include "HistMan/HistMan.h"
#include "AnalysisNtuples/ANtpDefaultValue.h"
                                                                                


CVSID("$Id: ValVtxModule.cxx,v 1.4 2008/07/11 19:19:48 boehm Exp $");

JOBMODULE(ValVtxModule, "ValVtxModule", "");
//......................................................................

ValVtxModule::ValVtxModule():
   counter(0),
   passcounter(0),
   passcutcounter(0),
   failcounter(0),
   writecounter(0)
{}

//......................................................................

ValVtxModule::~ValVtxModule()
{
}

//......................................................................

JobCResult ValVtxModule::Reco(MomNavigator* mom)
{
//======================================================================
// Reads in sue-style ntuples from mom, calculates a bunch of useful
// nue quantites
//======================================================================
   MSG("ValVtxModule",Msg::kDebug)<<"In ValVtxModule::Reco"<<endl;

   if(counter%1000==0){      
      MSG("ValVtxModule",Msg::kInfo)<<"On entry "<<counter<<endl;
   }
   counter++;

   TObject *obj=0;

   TIter objiter = mom->FragmentIter();
   while((obj=objiter.Next())){
      VtxRecord *vr = static_cast<VtxRecord *>(obj);
      if(vr){
         MSG("ValVtxModule",Msg::kDebug)<<"Found a VtxRecord in MOM"<<endl;
      }
      else{
         MSG("ValVtxModule",Msg::kError)<<"Didn't find a VtxRecord in MOM"<<endl;
         continue;
      }

      FillPlots(vr);
  }

  return JobCResult::kPassed; // kNoDecision, kFailed, etc.       
}

//......................................................................

void ValVtxModule::BeginJob()
{
  MSG("ValVtxModule",Msg::kDebug)<<"In ValVtxModule::BeginJob"<<endl;
                           
  names.clear();
  names.push_back("vfZError");   
  names.push_back("vfUError");
  names.push_back("vfVError");
  names.push_back("vfError");
  names.push_back("BestZError");
  names.push_back("BestUError");
  names.push_back("BestVError");
  names.push_back("BestError");
  names.push_back("relZError");
  names.push_back("relUError");
  names.push_back("relVError");
  names.push_back("relError");
                                                            
  classes.clear();
  classes.push_back("_far_nue");
  classes.push_back("_far_numu");
  classes.push_back("_far_bnue");
  classes.push_back("_far_nutau");
  classes.push_back("_far_nc");
  classes.push_back("_near_bnue");
  classes.push_back("_near_numu");
  classes.push_back("_near_nc");
  classes.push_back("_near_rock");
  classes.push_back("_far_rock");
  classes.push_back("_near_cosmic");
  classes.push_back("_far_cosmic");
  classes.push_back("_unknown");
                                                                                
  Int_t nbins = 1000;
  Float_t beg = -25.0;
  Float_t end = 25.0;                                              
  
  Int_t nbins2 = 250;
  Float_t beg2 = -5.0;
  Float_t end2 = 5.0;


                                
  for(UInt_t i = 0; i < names.size(); i++){
     for(UInt_t l = 0; l < classes.size(); l++){
           string dirstring = "vtxstats";
           HistMan hm2(dirstring.c_str());
           TString param = names[i]+classes[l];
//           hm2.Book<TH1F>(param,param,nbins,beg,end);
           if(i%4 == 0 || i%4 == 3) hm2.Book<TH1F>(param,param,nbins,beg,end);
           else hm2.Book<TH1F>(param,param,nbins,beg2,end2);
    }
  }

  MSG("ValVtxModule",Msg::kDebug)<<"Created "<<names.size()*classes.size()
           <<" histogram entries in BeginJob"<<endl;                
}

void ValVtxModule::EndJob()
{

   MSG("ValVtxModule",Msg::kInfo)<<"Counter "<<counter
                              <<" passcutcounter "<<passcutcounter
                              <<" passcounter "<<passcounter
                              <<endl;

   TFile* file = new TFile("Results.root", "update");
   file->cd();
   HistMan hm2("vtxstats");
   hm2.WriteOut(*file);

  delete file;

  MSG("ValVtxModule",Msg::kDebug)<<"File output complete"<<endl;
                                                                                       
}

void ValVtxModule::FillPlots(VtxRecord *vr)
{
     TString id = MakeIdString(vr);

     vector<Float_t> vals;
     vals.push_back(vr->vtxfind.vfErrZ);
     vals.push_back(vr->vtxfind.vfErrU);
     vals.push_back(vr->vtxfind.vfErrV);
     vals.push_back(vr->vtxfind.vfError);
     vals.push_back(vr->vtxfind.vfBestErrZ);
     vals.push_back(vr->vtxfind.vfBestErrU);
     vals.push_back(vr->vtxfind.vfBestErrV);
     vals.push_back(vr->vtxfind.vfBestError);
     vals.push_back(vr->vtxfind.vfRelErrZ);
     vals.push_back(vr->vtxfind.vfRelErrU);
     vals.push_back(vr->vtxfind.vfRelErrV);
     vals.push_back(vr->vtxfind.vfRelError);

     HistMan hm2("vtxstats");

     for(unsigned int i = 0; i < vals.size(); i++) {
        if(!ANtpDefVal::IsDefault(vals[i])){
             TString hname = names[i]+id;
             TH1F* hist = hm2.Get<TH1F>(hname);
             hist->Fill(vals[i], 1.0);
        }
     }
     MSG("ValVtxModule",Msg::kDebug)<<"Done filling Histograms"<<endl;
}

TString ValVtxModule::MakeIdString(VtxRecord *vr)
{
  Detector::Detector_t d = vr->GetHeader().GetVldContext().GetDetector();
                                                                                
  TString det, dm;
  TString type;
                                                                                
  if(d==Detector::kFar){
    det = "_far";
  }
  else if(d==Detector::kNear){
    det = "_near";
  }
  else{
    return "_unknown";
  }
                                                                                
  if(vr->mctrue.fNueClass == 0){type = "nc";}
  else if(vr->mctrue.fNueClass == 1){type = "numu";}
  else if(vr->mctrue.fNueClass == 2){type = "nue";}
  else if(vr->mctrue.fNueClass == 3){type = "nutau";}
  else if(vr->mctrue.fNueClass == 4){type = "bnue";}
  else { return "_unknown";}

  if(vr->mctrue.nuVtxZ < 0 && vr->mctrue.fNueClass == 1)  
       type = "rock";

  TString id = det+"_" + type;

  return id;
}

---