NtpTools.cxx

Go to the documentation of this file.

#include "NtpTools.h"

#include "BeamDataUtil/BeamMonSpill.h"
#include "BeamDataUtil/BMSpillAna.h"
#include "MessageService/MsgService.h"
#include "Conventions/Munits.h"
#include "EnergyCorrections.h"
#include <set>
#include <cassert>

#include "TruthHelperNtuple/NtpTHEvent.h"
#include "MCNtuple/NtpMCStdHep.h"

#include "MCReweight/MCReweight.h"
// #include "MCReweight/NeugenWeightCalculator.h"
// #include "MCReweight/ReweightHelpers.h"
#include "MCReweight/Zfluk.h"

// March 12 2008- (Pedro) Added security to avoid fpe when taking sqrt of very high CorrectSignedMomentumFromCurvature quantity. Happened only extremely rarely. 

CVSID("$Id: NtpTools.cxx,v 1.5 2009/05/27 01:01:09 rbpatter Exp $");

NtpTools::NtpTools()
  : fCanvas(NULL),
    fUZAxis(NULL),
    fVZAxis(NULL)
{
  MSG("NtpAna",Msg::kDebug) << " NtpTools::NtpTools() "<<endl;
  fCanvas = new TCanvas("canv", " NtpTools ", 800, 422);
  fFluxModel = new Zfluk();
  fFluxModel->UseParameterization("SKZP");  
}


NtpTools::~NtpTools(){
  if(fCanvas) delete fCanvas;
  if(fUZAxis) delete fUZAxis;
  if(fVZAxis) delete fVZAxis;
}

float NtpTools::GetMCEventWeight(const NtpStRecord* /*record*/, const int /*eventno*/){
  //
  //Calculates the SKZP weight for each event
  //

  /*
    const NtpMCTruth* mc = GetTruth(record, eventno);

    double par[]={1.5686,
    -2.9241,
    1.0259,
    1.1673,
    0.9719,
    2.4302,
    0.80088};
    float pz = mc->flux.tpz;
    float pt = sqrt(mc->flux.tpx*mc->flux.tpx+mc->flux.tpy*mc->flux.tpy);

    float weight = fFluxModel->GetWeight(mc->flux.tptype,pt,pz,par);
    return weight;
  */

  cout << "CAUTION: RETURNING 1.0 AS ZFLUK WEIGHT" << endl;

  return 1.0;

}

//----------------------------------------------------------------------------------
// Adapted from NueAna::StdHepInfoAna ! Gotta make sure definition hasn't changed ! 
// This variable is used for pt reweighting of hadron showers in nue analysis
//
// Pedro Ochoa
//
// May 11 2007
//--------------------------------------------------------------------------------- 
float NtpTools::GetTotHadPt(const NtpStRecord* record, const NtpMCTruth* truth){

  float totpt=0;
  int index=truth->index;

  vector<int> pid;
  vector<double> px;
  vector<double> py;
  vector<double> pz;
  vector<double> eng;
  vector<double> mass;

  double hspx = 0;
  double hspy = 0;
  double hspz = 0;
  double hse  = 0;

  //cout<<"counting pi0"<<endl;
  Int_t nStdHep = record->mchdr.nstdhep;
  NtpMCStdHep *ntpStdHep = 0;
  for(int i=0; i<nStdHep; i++){
    ntpStdHep = dynamic_cast<NtpMCStdHep *>(record->stdhep->At(i));
    //check that stdhep entry corresponds to this mc event
    if(ntpStdHep->mc!=index) continue;
    
    if(ntpStdHep->IstHEP==3){//resonance decay
      if (TMath::Abs(ntpStdHep->p4[3])>hse){
        hspx = ntpStdHep->p4[0];
        hspy = ntpStdHep->p4[1];
        hspz = ntpStdHep->p4[2];
        hse  = TMath::Abs(ntpStdHep->p4[3]);
      }
    }
    else {
      //check if this particle comes from a resonance decay or DIS
      bool resdecay = true;
      for (int ip = ntpStdHep->parent[0]; ip<=ntpStdHep->parent[1]; ip++){
        if (ip==-1) {
          resdecay = false;
          break;
        }
        NtpMCStdHep *parent = dynamic_cast<NtpMCStdHep *>(record->stdhep->At(ip));
        resdecay = resdecay && (parent->IstHEP==3);
        if (!resdecay) break;
      }
      if (resdecay&&ntpStdHep->IdHEP) {
        pid.push_back(ntpStdHep->IdHEP);
        px.push_back(ntpStdHep->p4[0]);
        py.push_back(ntpStdHep->p4[1]);
        pz.push_back(ntpStdHep->p4[2]);
        eng.push_back(ntpStdHep->p4[3]);
        mass.push_back(ntpStdHep->mass);
      }
    }
  }

  double Phs = sqrt(pow(hspx,2)+pow(hspy,2)+pow(hspz,2));
  for (unsigned ipar = 0; ipar<pid.size(); ipar++){//loop through hadrons
    double ptot = sqrt(pow(px[ipar],2)+pow(py[ipar],2)+pow(pz[ipar],2));
    double Pz = (px[ipar]*hspx+py[ipar]*hspy+pz[ipar]*hspz)/Phs;
    double pt = pow(ptot,2)-pow(Pz,2);
    if (pt>0.00000001) pt = sqrt(pt);
    else pt = 0;

    totpt += pt;
  }

  return totpt;

}

bool NtpTools::FillRecoEInfo( const NtpStRecord* record, const int eventno, RecoE* e_reco){
  e_reco->Reset();
  const NtpSRTrack* track = GetPrimaryTrack(record, eventno);
  const NtpSREvent* event = GetEvent(record, eventno);
  Detector::Detector_t detector =record->GetHeader().GetVldContext().GetDetector();
  const     Bool_t isdata= record->GetHeader().GetVldContext().GetSimFlag()==SimFlag::kData;  
  const int pitt_flag = GetPittContainmentFlag(event);
  e_reco->EOkay = 1;

  if(track){
    if(pitt_flag ==1 || pitt_flag==3){
      e_reco->trkEMethod = 1;
    }else if(pitt_flag==2 || pitt_flag==4){
      e_reco->trkEMethod = 2;
    }
    Float_t mr=CorrectMomentumFromRange(track->momentum.range,isdata, detector);
    e_reco->trkERange = sqrt(mr*mr+0.105*0.105);
  
    if (track->momentum.qp!=0){
      Float_t p= CorrectSignedMomentumFromCurvature(1./track->momentum.qp,isdata,detector);
      if(p>0 && p<10000){
        e_reco->trkECurve =  sqrt(pow(p,2)+pow(0.105,2));
      } else e_reco->trkECurve = 0;
    }
    else{
      e_reco->trkECurve = 0;
      e_reco->EOkay = 0;
    }
    if(e_reco->trkEMethod==1) e_reco->trkE = e_reco->trkERange;
    else if(e_reco->trkEMethod==2)  e_reco->trkE = e_reco->trkECurve;
    else e_reco->EOkay = 0;
  }else{
    e_reco->trkE = 0;
    e_reco->trkEMethod = 0;
    e_reco->trkERange = 0;
    e_reco->trkECurve = 0;   
  }//tracks

  //now add in the other showers
  e_reco->shwE = 0;
  for(int ishw = 0; ishw<event->nshower; ++ishw){
    const NtpSRShower* ishower = dynamic_cast<const NtpSRShower*>(record->shw->At(event->shw[ishw]));
    float showere = CorrectShowerEnergy(ishower->shwph.linCCgev, 
                                        detector, 
                                        CandShowerHandle::kCC, 
                                        1, isdata);//or kWtCC but don't use that
    assert(ishower);
    e_reco->shwE+=showere;
  }
  
  e_reco->totalE= (e_reco->trkE + e_reco->shwE);
  return true;
}

int NtpTools::GetNEventsPerSlice(const NtpStRecord* record, const int eventno){
  const NtpSREvent* event = GetEvent(record, eventno);
  if(!event) return -1;

  int nmatch = 0;
  for(int ievent = 0; ievent<record->evthdr.nevent; ++ievent){
    const NtpSREvent* tmpevent = dynamic_cast<const NtpSREvent*>(record->evt->At(ievent));
    if(tmpevent->slc == event->slc){
      nmatch++;
    }
  }
  return nmatch;
}

Int_t NtpTools::GetPittContainmentFlag(const NtpSREvent* evt)
{
 //Calc Pittsburgh (Donna/Debdatta) containment flag

  Int_t pitt_flag=0;
  // pitt_flag:
  // 1 = track is fully     contained in the upstream   region
  // 2 = track is partially contained in the upstream   region
  // 3 = track is fully     contained in the downstream region
  // 4 = track is partially contained in the downstream region

  Float_t x=evt->end.x;
  Float_t y=evt->end.y;
  Float_t z=evt->end.z;
  Float_t u=evt->end.u;
  Float_t v=evt->end.v;
  Float_t rad=TMath::Sqrt(x*x + y*y);
  
  Float_t calZ=7;
  Float_t specZ=15.6;//same as Mad
  Float_t minR=0.8;

  Bool_t down_stop=false;
  Bool_t down_exit=false;
  Bool_t up_stop=false;
  Bool_t up_exit=false;
  
  // downstream
  if (z>calZ) {
    if (u>-0.85 && u<2.1 && 
        v>-2.1 && v<0.85 && 
        x>-1.5 && x<2.4 && 
        y>-1.4 && y<1.4  && 
        rad>minR && 
        z>calZ && z<=specZ) down_stop=true;
    else down_exit=true;
  }
  //ustream
  else if (z<=calZ && z>=0) {
    if (u>0.3 && u<1.8 && 
        v>-1.8 && v<-0.3 && 
        x<2.4 && 
        rad>minR &&
        z<calZ && z>=0) up_stop=true;
    else up_exit=true;
  }
  else cout<<"Ahhh, z="<<z<<endl;

  // set the track's pitt flag
  if (down_stop) pitt_flag = 3;
  if (up_stop  ) pitt_flag = 1;
  if (down_exit) pitt_flag = 4;
  if (up_exit  ) pitt_flag = 2;

  return pitt_flag;
}

bool NtpTools::FillMyShwInfo(const NtpStRecord* record, const int eventno, MyShwInfo* shwinfo){
  shwinfo->Reset();
  static float planePH[500];
  for(int i=0; i<500; ++i) planePH[i] = 0;
  const NtpSREvent* event = GetEvent(record, eventno);
  for(int istp = 0; istp<event->nstrip; ++istp){
    const NtpSRStrip* strip = dynamic_cast<const NtpSRStrip*>(record->stp->At(event->stp[istp]));
    shwinfo->totalPH+=(strip->ph0.sigcor+strip->ph1.sigcor);
    if(strip->plane<500){
      planePH[strip->plane]+=(strip->ph0.sigcor+strip->ph1.sigcor);
    }else{
      cout << "Failed to get plane number; " << strip->plane<< endl;
    }
  }
  shwinfo->maxPH3Plane = 0;
  for(int i=1; i<499; ++i){
    float tmp = planePH[i-1] +  planePH[i+1] +  planePH[i];
    if(tmp>shwinfo->maxPH3Plane) shwinfo->maxPH3Plane = tmp;
  }
  return true;
}

const NtpSRTrack* NtpTools::GetPrimaryTrack(const NtpStRecord* record, const int eventno){
  const NtpSREvent* event = dynamic_cast<const NtpSREvent*>(record->evt->At(eventno));
  if(event==NULL){
    MSG("NtpAna",Msg::kError) << " NtpTools::GetPrimaryTrack passed null event " <<endl;
  }
  int best_track = -1;
  int nplanes = 0;
  for(int itrk=0; itrk<event->ntrack; ++itrk){
    const NtpSRTrack* track = dynamic_cast<const NtpSRTrack*>(record->trk->At(event->trk[itrk]));

    int tpln = abs(track->end.plane - track->vtx.plane);
    if(tpln>nplanes) {
      nplanes = tpln;
      best_track = itrk;
    }
  }
  if(best_track==-1){
    return NULL;
  }else{
    return (dynamic_cast<const NtpSRTrack*>(record->trk->At(event->trk[best_track])));
  }
}

bool NtpTools::PassBeamCuts(const BeamType_t beamtype, const NtpBDLiteRecord* beamntp, float& pot, bool isMC ){
  //
  //PassBeamCuts
  //

  if(isMC){
   pot =  24.2;
   return true;
  }

  static Registry* beam_config = NULL;
  static BeamMonSpill spill;
  static BMSpillAna bmsa;
  static BeamType_t lastbeamtype = beamtype;

  if(beamntp==NULL) return 0;
  
  if(beam_config==NULL){
    beam_config = new Registry();
    beam_config->UnLockValues();
    //beam_config->Set("TargetIn",1);
    //beam_config->Set("BeamType",-1);//don't care if ME,LE,etc
    //beam_config->Set("PosTgtXMin",-2.0*Munits::mm);
    beam_config->Set("BeamType",(int)(beamtype));
    //have to make sure that zero is not included
    beam_config->Set("PosTgtXMax",-0.0000001*Munits::mm);
    beam_config->Set("PosTgtYMin", 0.0000001*Munits::mm);
    //beam_config->Set("PosTgtYMax", 2.0*Munits::mm);
    beam_config->LockValues();
    bmsa.Config(*beam_config);    
    //print the registry
    beam_config->PrettyPrint(cout);
  }
  if(lastbeamtype!=beamtype){
    beam_config->UnLockValues();
    beam_config->Set("BeamType",(int)(beamtype));
    beam_config->LockValues();
    bmsa.Config(*beam_config);    
    lastbeamtype = beamtype;
  }
  
  bmsa.SetSpill(*beamntp,spill);
  pot =  beamntp->tortgt;
  if (!bmsa.SelectSpill()){
    MAXMSG("NtpAna",Msg::kInfo,10)
      <<"Rejected spill!"
      <<", tortgt="<<beamntp->tortgt
      <<", timeDiff="<<beamntp->GetHeader().GetTimeDiffStreamSpill() <<endl;
    return 0;    
  }else{
    return true;
  }
}

//Drawing stuff

bool NtpTools::ResetCanvas(Detector::Detector_t det){
  if(fUZAxis) delete fUZAxis;
  if(fVZAxis) delete fVZAxis;
  fUZAxis = NULL;
  fVZAxis = NULL;
  fUZMarkers.clear();
  fVZMarkers.clear();
  fCanvas->Clear();
  fCanvas->Divide(1,2);
  
  switch(det){
  case Detector::kFar:
    fUZAxis = new TH2F("fUZAxis", "", 100, 0, 17, 100, -4, 4.);
    fVZAxis = new TH2F("fVZAxis", "", 100, 0, 17, 100, -4, 4.);    
    break;
  case Detector::kNear:
    fUZAxis = new TH2F("fUZAxis", "", 100, 0, 17, 100, -1.5, 3.);
    fVZAxis = new TH2F("fVZAxis", "", 100, 0, 17, 100, -3, 1.5);    
    break;
  default:
    return false;
    break;
  };
  
  fUZAxis->SetXTitle("Z/m");
  fUZAxis->SetYTitle("U/m");
  fVZAxis->SetXTitle("Z/m");
  fVZAxis->SetYTitle("V/m");
  return true;
}

bool  NtpTools::MakeEventPict(const char* title, const NtpStRecord* record, const int eventno){
  bool rc = false;
  if(ResetCanvas(GetDetector(record))==false) return false;
  
  const NtpSREvent* event = GetEvent(record, eventno);
  const TClonesArray* strips = (record->stp);
  const TClonesArray* tracks = (record->trk);
  
  std::set<int> trkstp_index;
                                                
 for(int itrk=0; itrk<event->ntrack; ++itrk){
    const NtpSRTrack* track = dynamic_cast<const NtpSRTrack*>(tracks->At(event->trk[itrk]));
    if(track==NULL){
      MSG("NtpAna",Msg::kError) << " failed to cast track " <<endl;
    }else{
      for(int istp=0; istp<track->nstrip; istp++){
        trkstp_index.insert(track->stp[istp])      ;
      }    
    }
  }

  for(int istp =0; istp< event->nstrip ; ++istp){
    const int stp_index = event->stp[istp];
    const NtpSRStrip* strip = dynamic_cast<const NtpSRStrip*>(strips->At(stp_index));
    if(strip==NULL){
      MSG("NtpAna",Msg::kError) << " failed to cast strip " <<endl;
      rc = false;
    }else{
      
      TMarker tmpmarker(strip->z, strip->tpos, 
                        20);
      tmpmarker.SetMarkerSize(0.5);
      if(trkstp_index.count(stp_index)>0)      tmpmarker.SetMarkerColor(2);
      if(strip->planeview==2){
        fUZMarkers.push_back(tmpmarker);
      }else{
        fVZMarkers.push_back(tmpmarker);
      }
    }
    
  }

  // RBP - changed "!=" to "=!" in a guess at what was intended by the original author.
  rc = !PrintCanvas(title);
  
  return rc;
}

bool NtpTools::PrintCanvas(const char* title){
  fCanvas->cd(1);
  fUZAxis->Draw("axis");
  float minz = 999;
  float maxz = 0;

  for(unsigned int i=0; i<fUZMarkers.size(); i++){
    fUZMarkers[i].Draw();
    if(minz>fUZMarkers[i].GetX()) minz=fUZMarkers[i].GetX();
    if(maxz<fUZMarkers[i].GetX()) maxz=fUZMarkers[i].GetX();
  }
  fUZAxis->GetXaxis()->SetRangeUser(minz, maxz);

  fCanvas->cd(2);
  fVZAxis->Draw("axis");
  minz = 999;
  maxz = 0;

  for(unsigned int i=0; i<fVZMarkers.size(); i++){
    fVZMarkers[i].Draw();
    if(minz>fVZMarkers[i].GetX()) minz=fVZMarkers[i].GetX();
    if(maxz<fVZMarkers[i].GetX()) maxz=fVZMarkers[i].GetX();
  }
  fVZAxis->GetXaxis()->SetRangeUser(minz, maxz);
  
  fCanvas->Print(title);
  return true;
}

const NtpMCTruth* NtpTools::GetTruth(const NtpStRecord* record, const int eventno){
    if(!record || record->thevt==NULL){
    cout << "*** no truth helper! " <<endl;
    return NULL;
  }
  if(record->evt->GetEntries() != record->thevt->GetEntries()){
    cout << " funny matching evt - >thevt " <<endl;
    return NULL;
  }
  const NtpTHEvent* thevent = dynamic_cast<const NtpTHEvent*>(record->thevt->At(eventno));
  if(!thevent){
    cout << "*** unable to get thevt" <<endl;
    return NULL;
  }
  return (dynamic_cast<const NtpMCTruth*>(record->mc->At(thevent->neumc)));
  
}

bool NtpTools::GetTrueEnergy(const NtpStRecord* record, const int eventno, float& mceneu){
  const NtpMCTruth* mcevent = GetTruth(record, eventno);

  if(!mcevent){
    cout << "*** unable to get mcevent" <<endl;
    return false;    
  }
  mceneu = fabs(mcevent->p4neu[3]);
  return true;
}

---