ANtpInfoObjectFillerBeam.cxx

Go to the documentation of this file.

//$Id: ANtpInfoObjectFillerBeam.cxx,v 1.40 2010/01/19 18:10:24 tinti Exp $
//
//ANtpInfoObjectFiller.cxx
//
//B. Rebel 02/2005

#include "AnalysisNtuples/Module/ANtpInfoObjectFillerBeam.h"
#include "MessageService/MsgService.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "StandardNtuple/NtpStRecord.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "BField/BfieldCoilCurrent.h"
#include "BeamDataUtil/BDSpillAccessor.h"
#include "BeamDataUtil/BeamMonSpill.h"
#include "BeamDataNtuple/NtpBDLiteRecord.h"
#include "SpillTiming/SpillTimeFinder.h"
#include "SpillTiming/SpillServerMon.h"
#include "SpillTiming/SpillServerMonFinder.h"
#include "TStopwatch.h"
#include "Conventions/BeamType.h"
#include <cassert>
#include <algorithm>

using namespace BeamType;

ClassImp(ANtpInfoObjectFillerBeam)

CVSID("$Id: ANtpInfoObjectFillerBeam.cxx,v 1.40 2010/01/19 18:10:24 tinti Exp $");

//-------------------------------------------------------------------
ANtpInfoObjectFillerBeam::ANtpInfoObjectFillerBeam()
{
  
  fSpillAna = new BMSpillAna();
  //use the database to determine the cuts for this data sample
  fSpillAna->UseDatabaseCuts();
  
  MSG("ANtpInfoObjectFillerBeam", Msg::kDebug) << "ANtpInfoObjectFillerBeam::Constructor" << endl;
  
}

//-------------------------------------------------------------------
ANtpInfoObjectFillerBeam::ANtpInfoObjectFillerBeam(Detector::Detector_t detector) :
  ANtpInfoObjectFiller(detector)
{

  fSpillAna = new BMSpillAna();
  //use the database to determine the cuts for this data sample
  fSpillAna->UseDatabaseCuts();
    
  MSG("ANtpInfoObjectFillerBeam", Msg::kDebug) << "ANtpInfoObjectFillerBeam::Constructor" << endl;
   
}

//----------------------------------------------------------------------
ANtpInfoObjectFillerBeam::~ANtpInfoObjectFillerBeam()
{

  delete fSpillAna;
  MSG("ANtpInfoObjectFillerBeam", Msg::kDebug) << "ANtpInfoObjectFillerBeam::Destructor" << endl;
  
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to beam monitoring information
void ANtpInfoObjectFillerBeam::FillBeamInformation(VldTimeStamp &timeStamp,
                                                   Detector::Detector_t det,
                                                   SimFlag::SimFlag_t dataType,
                                                   ANtpBeamInfo *beamInfo)
{
  //get the appropriate context
  VldContext vldc(det, dataType, timeStamp);

  beamInfo->Reset();

//   TStopwatch stopwatch;
//     stopwatch.Start();
//     stopwatch.Stop();
//     MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime()
//                                            << " to make LoadSpill" << endl;
//     stopwatch.Reset();

  if(dataType == SimFlag::kData){
    //make the thing to access the beam information from the db
    BDSpillAccessor& spillAccess = BDSpillAccessor::Get();

    double x = 0.;
    double y = 0.;
    double xrms = 0.;
    double yrms = 0.;

    //load the spill for the validity context
    const BeamMonSpill &spillInfo = *(spillAccess.LoadSpill(timeStamp));
    fSpillAna->SetSpill(spillInfo);
    spillInfo.BpmAtTarget(x,y,xrms,yrms);

    beamInfo->tor101 = spillInfo.fTor101;
    beamInfo->tr101d = spillInfo.fTr101d;
    beamInfo->tortgt = spillInfo.fTortgt;
    beamInfo->trtgtd = spillInfo.fTrtgtd;
    beamInfo->hornCurrent = spillInfo.fHornCur;
    beamInfo->targetBPMX = x;
    beamInfo->targetBPMY = y;
    beamInfo->profileWidthX = spillInfo.fProfWidX;
    beamInfo->profileWidthY = spillInfo.fProfWidY;
    beamInfo->hadronIntensity = spillInfo.fHadInt;
    beamInfo->muonIntensity1 = spillInfo.fMuInt1;
    beamInfo->muonIntensity2 = spillInfo.fMuInt2;
    beamInfo->muonIntensity3 = spillInfo.fMuInt3;
    beamInfo->beamType = (Int_t)( spillInfo.BeamType() );
    beamInfo->protonIntensity = beamInfo->tortgt;

    //cout << BeamType::AsString(spillInfo.BeamType()) << " " << beamInfo->beamType << endl;
    
    SpillTimeFinder &stf = SpillTimeFinder::Instance();
    VldTimeStamp tons = stf.GetTimeOfNearestSpill(vldc);
    beamInfo->timeToNearestSpill = stf.GetTimeToNearestSpill(vldc);
    beamInfo->nearestSecToSpill = tons.GetSec();
    beamInfo->nearestNSToSpill = tons.GetNanoSec();
    VldTimeStamp ds = spillInfo.SpillTime()-vldc.GetTimeStamp();
    fSpillAna->SetTimeDiff(ds.GetSeconds());
    beamInfo->streamSpillTimeDiff = ds.GetSeconds();
    beamInfo->goodSpill = (int)fSpillAna->SelectSpill();

    //gps error is only relevant for the far detector, so 
    //only ping the database for the far detector
    if(det == Detector::kFar){
      SpillServerMonFinder& spillFinder = SpillServerMonFinder::Instance();
      const SpillServerMon& nearestSpill= spillFinder.GetNearestSpill(vldc);
      tons = nearestSpill.GetSpillTime()-timeStamp;
      beamInfo->deltaSecToSpillGPS = TMath::Abs(tons.GetSec());
      beamInfo->gpsError = nearestSpill.GetSpillTimeError();
    }
  }
  else if(dataType == SimFlag::kMC){
//     beamInfo->tor101 = spillInfo->fTor101;
//     beamInfo->tortgt = spillInfo->fTortgt;
//     beamInfo->trtgtd = spillInfo->fTrtgtd;
//     beamInfo->hornCurrent = spillInfo->fHornCur;
//     beamInfo->targetBPMX = spillInfo->fTargBpmX[0];
//     beamInfo->targetBPMY = spillInfo->fTargBpmY[0];
//     beamInfo->profileWidthX = spillInfo->fProfWidX;
//     beamInfo->profileWidthY = spillInfo->fProfWidY;
//     beamInfo->hadronIntensity = spillInfo->fHadInt;
//     beamInfo->muonIntensity1 = spillInfo->fMuInt1;
//     beamInfo->muonIntensity2 = spillInfo->fMuInt2;
//     beamInfo->muonIntensity3 = spillInfo->fMuInt3;
//     beamInfo->targetPosition = (int)spillInfo->GetStatusBits().beam_type;
//     beamInfo->protonIntensity = beamInfo->tor101;
  }
  
  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to beam monitoring information
void ANtpInfoObjectFillerBeam::FillBeamInformation(NtpBDLiteRecord *record,
                                                   ANtpBeamInfo *beamInfo)
{

  beamInfo->Reset();

  //got to get the header to get the status bits
  if(record){

    double x = 0.;
    double y = 0.;
    double xrms = 0.;
    double yrms = 0.;

    //need a BeamMonSpill object to use the BMSpillAna object
    BeamMonSpill spillInfo;

    //set the spill in the BMSpillAna object
    fSpillAna->SetSpill(*record, spillInfo);
    spillInfo.BpmAtTarget(x,y,xrms,yrms);

    beamInfo->tor101 = spillInfo.fTor101;
    beamInfo->tr101d = spillInfo.fTr101d;
    beamInfo->tortgt = spillInfo.fTortgt;
    beamInfo->trtgtd = spillInfo.fTrtgtd;
    beamInfo->hornCurrent = spillInfo.fHornCur;
    beamInfo->targetBPMX = x;
    beamInfo->targetBPMY = y;
    beamInfo->profileWidthX = spillInfo.fProfWidX;
    beamInfo->profileWidthY = spillInfo.fProfWidY;
    beamInfo->hadronIntensity = spillInfo.fHadInt;
    beamInfo->muonIntensity1 = spillInfo.fMuInt1;
    beamInfo->muonIntensity2 = spillInfo.fMuInt2;
    beamInfo->muonIntensity3 = spillInfo.fMuInt3;
    beamInfo->beamType = (Int_t)( spillInfo.BeamType() );
    beamInfo->protonIntensity = beamInfo->trtgtd;
    
    beamInfo->timeToNearestSpill = record->dt_nearest;
    beamInfo->streamSpillTimeDiff = record->GetHeader().GetTimeDiffStreamSpill();
    beamInfo->nearestSecToSpill = record->nearest_sec;
    beamInfo->nearestNSToSpill = record->nearest_nsec;
    beamInfo->goodSpill = (int)fSpillAna->SelectSpill();
  }
  
  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to the Beam MC
void ANtpInfoObjectFillerBeam::FillBeamMCTruthInformation(NtpMCTruth *ntpMCTruth, 
                                                          TClonesArray *stdArray,
                                                          ANtpTruthInfoBeam *truthInfo)
{  
  //the NtpMCTruth object stdhep array tells you which entries in the stdhep array 
  //are the first and last for the event
  Int_t lowIndex = (Int_t)ntpMCTruth->stdhep[0];
  Int_t highIndex = (Int_t)ntpMCTruth->stdhep[1];

  MSG("ANtpInfoObjectFillerBeam", Msg::kDebug) << "fill beam mc truth " 
                                              << lowIndex << " " << highIndex
                                              << endl;

  //already fill this variable in the FillMCTruthInformation() method
  //but leave it in to quiet the compiler
  truthInfo->nuFlavor = ntpMCTruth->inu;
  truthInfo->nonOscNuEnergy = ntpMCTruth->p4neunoosc[3];
  truthInfo->nonOscNuDCosX = ntpMCTruth->p4neunoosc[0];
  truthInfo->nonOscNuDCosY = ntpMCTruth->p4neunoosc[1];
  truthInfo->nonOscNuDCosZ = ntpMCTruth->p4neunoosc[2];
  truthInfo->nonOscNuFlavor = ntpMCTruth->inunoosc;
  truthInfo->resonanceCode = ntpMCTruth->iresonance;
  
  truthInfo->atomicWeight = ntpMCTruth->a;
  truthInfo->atomicNumber = ntpMCTruth->z;
  truthInfo->bjorkenX = ntpMCTruth->x;
  truthInfo->q2 = ntpMCTruth->q2;
  truthInfo->w2 = ntpMCTruth->w2;
  truthInfo->sigma = ntpMCTruth->sigma;
  truthInfo->emShowerFraction = ntpMCTruth->emfrac;

  truthInfo->targetExitX = ntpMCTruth->flux.tvx;
  truthInfo->targetExitY = ntpMCTruth->flux.tvy;
  truthInfo->targetExitZ = ntpMCTruth->flux.tvz;
  truthInfo->targetParentPX = ntpMCTruth->flux.tpx;
  truthInfo->targetParentPY = ntpMCTruth->flux.tpy;
  truthInfo->targetParentPZ = ntpMCTruth->flux.tpz;
  truthInfo->targetParentType = ntpMCTruth->flux.tptype;
  truthInfo->parentX = ntpMCTruth->flux.vx;
  truthInfo->parentY = ntpMCTruth->flux.vy;
  truthInfo->parentZ = ntpMCTruth->flux.vz;
  truthInfo->parentPX = ntpMCTruth->flux.pdpx;
  truthInfo->parentPY = ntpMCTruth->flux.pdpy;
  truthInfo->parentPZ = ntpMCTruth->flux.pdpz;
  truthInfo->parentPID = ntpMCTruth->flux.ptype;

  NtpMCStdHep *ntpMCStdHep = 0;
  
  //loop over the entries to figure out the targert nucleon
  Int_t nucleon = 0;
  Int_t nuFlavor = ntpMCTruth->inu;

  //gotta add up Mike K.'s true visible energy - email from mike on 5/26/06
  //Brian: You must define "true visible energy" for our regular MC files.  
  //It is not y*Enu for you.  You do it taking the *kinetic energy* of all
  //baryons and the *total energy* (e.g. include mass) of mesons. Use only
  //particles with IstHEP==1. Do add in any photons (e.g., from nucleon
  //resonance decays).  Use the code on all events, even those with zero
  //shower energy. Let me know if you have questions or problems.
  //in practice for everything that has 16 < abs(IdHep) < 1e5 and is not
  //2212 (proton) or 2112 (neutron) add total energy (p4[3]).  if is a 
  //proton or neutrino add sqrt(p4[3]*p4[3] - mass*mass)
  
  truthInfo->trueVisibleE = 0.;

  int idHEP = 0;
  int istHEP = 0;
  bool gotOneAlreadyMate = false;  //used for HadronicFinalState

  for(Int_t i = lowIndex; i < highIndex+1; ++i){
    
    ntpMCStdHep = dynamic_cast<NtpMCStdHep *>(stdArray->At(i));

    idHEP = ntpMCStdHep->IdHEP;
    istHEP = ntpMCStdHep->IstHEP;

    if(TMath::Abs(idHEP) > 16 
       && TMath::Abs(idHEP) < 1e5 
       && istHEP == 1){
      
      if(idHEP == 2212 || idHEP == 2112) 
        truthInfo->trueVisibleE += ntpMCStdHep->p4[3] - ntpMCStdHep->mass;
      else
        truthInfo->trueVisibleE += ntpMCStdHep->p4[3];
    }
    
    if(istHEP == 11){
      if(idHEP == 2212) nucleon = 1;
      else if(idHEP == 2112) nucleon = 0;
      else if(TMath::Abs(idHEP) > 1000000000)  nucleon = 2;
      else nucleon = 3;
    }
    
    if(nucleon == 1 && nuFlavor > 0) truthInfo->initialState = 1;       //p + v
    else if(nucleon == 0 && nuFlavor > 0) truthInfo->initialState = 2;  //n + v
    else if(nucleon == 1 && nuFlavor < 0) truthInfo->initialState = 3;  //p + vbar
    else if(nucleon == 0 && nuFlavor < 0) truthInfo->initialState = 4;  //n + vbar
    else if(nucleon == 2 && nuFlavor > 0) truthInfo->initialState = 5;  //N + v
    else if(nucleon == 3 && nuFlavor > 0) truthInfo->initialState = 6;  //A + v
    else if(nucleon == 2 && nuFlavor < 0) truthInfo->initialState = 7;  //N + vbar
    else if(nucleon == 3 && nuFlavor < 0) truthInfo->initialState = 8;  //A + vbar      

    //figure out the hadronic final state too, but only for non RES processes
    if(istHEP == 3 && !gotOneAlreadyMate){
      if(ntpMCTruth->iresonance != 1002 || TMath::Abs(idHEP) != 15){
        truthInfo->hadronicFinalState = idHEP%1000;
        gotOneAlreadyMate = true;
      }         
    }
    
  }//end loop over stdhep array for this event
  
  return;
}

---