NtpSRBleachFiller.cxx

Go to the documentation of this file.

//
// NtpSRBleachFiller
//
// Package: CandNtupleSR/Module
//
// Purpose: Methods for filling an NtpSRBleach object.
//
#include <iostream>
using namespace std;

#include "CandNtupleSR/Module/NtpSRBleachFiller.h"
#include "CandNtupleSR/Module/NtpSRNDAPPlaneHistory.h"
#include "CandData/CandRecord.h"
#include "RawData/RawRecord.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawDigit.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandEventHandle.h"
#include "Record/RecMinosHdr.h"
#include "MessageService/MsgService.h"
#include "Validity/VldContext.h"

CVSID("$Id: NtpSRBleachFiller.cxx,v 1.2 2005/11/07 03:04:02 schubert Exp $");

// Definition of methods (alphabetical order)
// ***************************************************

//......................................................................
double NtpSRBleachFiller::GetEventDuration(const CandEventHandle* cndevt) {
  //
  // Purpose: Calculate event duration from first and last strip time in 
  //          event. 
  //
  // Arguments: CandEventHandle ptr.
  // 
  // Return: event duration in nsec.
  //
  // Author: Tobias Raufer (tobias.raufer@physics.ox.ac.uk) 
  //
  // Creation Date: 11/05

  Double_t result = -1;
  if (!cndevt){
    MSG("NtpSR",Msg::kWarning) 
      << "No CandEventHandle Found, will return the default value"<<endl;
    return result;
  }

  const VldContext* vldc = cndevt->GetVldContext();
  Detector::Detector_t fDetType = vldc->GetDetector();
  if (fDetType!=Detector::kNear) return result;

  // loop over all strips and find earliest and latest strip time
  Double_t evtStart = 99999;
  Double_t evtEnd = -1;
  CandStripHandleItr stpItr(cndevt->GetDaughterIterator());
  while (CandStripHandle* strip=dynamic_cast<CandStripHandle*>(stpItr())){
   
    if (strip->GetTime() < evtStart) {
      evtStart = strip->GetTime();
    }
    if (strip->GetTime() > evtEnd) {
      evtEnd = strip->GetTime();
    }
 }
  result = (evtEnd - evtStart)*1e9; // result is in nsec
  return result;
}

//......................................................................
double NtpSRBleachFiller::GetFixedWindowPH(const CandEventHandle* cndevt, 
                                           const CandRecord* cndrec) {
  //
  // Purpose: Calculate predicted afterpulsing ph in a fixed window around
  //          the event vertex.
  //
  // Arguments: CandEventHandle ptr, CandRecord ptr.
  // 
  // Return: fixed window ph.
  //
  // Author: Tobias Raufer (tobias.raufer@physics.ox.ac.uk) 
  //
  // Creation Date: 11/05

  Double_t result = -1;
  if (!cndevt){
    MSG("NtpSR",Msg::kWarning) 
      << "No CandEventHandle Found, will return the default value"<<endl;
    return result;
  }
  if (!cndrec){
    MSG("NtpSR",Msg::kWarning) 
      << "No CandRecord Found, will return the default value"<<endl;
    return result;
  }

  VldContext vldc = cndrec->GetHeader()->GetVldContext();
  Detector::Detector_t fDetType = vldc.GetDetector();
  if (fDetType!=Detector::kNear) return result;

  // construct a NtpSRNDAPPlaneHistory object for each plane
  // this class contains the exponential model for the afterpulsing
  // based on the previous activity in the snarl
  map<Int_t,NtpSRNDAPPlaneHistory> planeHistMap; 
  
  for (Int_t plane=1; plane<282; plane++) {
    NtpSRNDAPPlaneHistory h(plane);
    planeHistMap.insert(std::pair<Int_t,NtpSRNDAPPlaneHistory>(plane,h));
  }

  // loop over all strips and fill NtpSRNDAPPlaneHistory objects
  const CandStripListHandle *striplisthandle 
    = dynamic_cast <const CandStripListHandle*> 
    (cndrec->FindCandHandle("CandStripListHandle"));
  if ( !striplisthandle ) {
    MSG("NtpSR",Msg::kWarning) 
      << "No strip list found in CandRecord. "
      << "Will return the default value " << endl;
    return result; // no strips => done
  }

  TIter stripItr(striplisthandle->GetDaughterIterator());
  while ( CandStripHandle* strip=dynamic_cast<CandStripHandle*>(stripItr())) {
    
    Int_t plane = strip->GetPlane();

    // add strip to NtpSRNDAPPlaneHistory object
    map<Int_t,NtpSRNDAPPlaneHistory>::iterator pos;    
    pos = planeHistMap.find(plane);
    if (pos!=planeHistMap.end())
      pos->second.AddStrip(strip->GetTime(),
                           strip->GetCharge(CalDigitType::kSigCorr));  
    else 
      MSG("NtpSR",Msg::kWarning) 
        << "No NtpSRNDAPPlaneHistory object for plane " << plane 
        << ". Not using this strip!" << endl;
  }

  // calculate pred. afterpulsing in 33 plane window around the vertex
  
  // find vertex plane; if track exists, use track vertex
  Int_t vtxPlane = -1; 
  const CandTrackHandle* trk = cndevt->GetPrimaryTrack();
  if (trk) 
    vtxPlane = trk->GetVtxPlane();
  else vtxPlane = cndevt->GetVtxPlane();
  
  Int_t minPlane = max(vtxPlane-2,1);
  Int_t maxPlane = min(vtxPlane+30,281);

  Double_t evtTime = 9999; // get event time from first strip in event
  CandStripHandleItr stpItr(cndevt->GetDaughterIterator());
  while (CandStripHandle* strip=dynamic_cast<CandStripHandle*>(stpItr())){
    if (strip->GetTime() < evtTime) {
      evtTime = strip->GetTime();
    }
  }

  Double_t APsum= 0;
  for (Int_t pl = minPlane; pl<maxPlane;pl++) {
    map<Int_t,NtpSRNDAPPlaneHistory>::iterator pos;    
    pos = planeHistMap.find(pl);
    if (pos!=planeHistMap.end())
      APsum +=
        pos->second.APPrediction(evtTime);
    else 
      MSG("NtpSR",Msg::kWarning) 
        << "No NtpSRNDAPPlaneHistory object for plane " << pl 
        << ". Not adding afterpulsing prediction for this plane!" 
        << endl;
  }
  result = APsum;
  return result;
}

//......................................................................
float NtpSRBleachFiller::GetlateBucketPHFraction(const CandEventHandle* cndevt,
                                                 const RawRecord* rawrec) {
  //
  // Purpose: Calculate late bucket ph fraction for each event.
  //
  // Arguments: CandEventHandle ptr, RawRecord ptr.
  // 
  // Return: late bucket ph fraction. -1 if error.
  //
  // Authors: Tingjun Yang (tjyang@stanford.edu) 
  //          Jiajie Ling  (ling@barney.physics.sc.edu)
  //
  // Creation date: 11/05 

  float ratio = -1;
  if (!cndevt){
    MSG("NtpSR",Msg::kWarning) 
      << "No CandEventHandle Found, will return the default value"<<endl;
    return ratio;
  }
  if (!rawrec){
    MSG("NtpSR",Msg::kWarning) 
      << "No RawRecord Found, will return the default value"<<endl;
    return ratio;
  }

  VldContext vldc = rawrec->GetHeader()->GetVldContext();
  Detector::Detector_t fDetType = vldc.GetDetector();
  if (fDetType!=Detector::kNear) return ratio;

  const RawDigitDataBlock *rddb = dynamic_cast<const RawDigitDataBlock *>
                                  (rawrec->FindRawBlock("RawDigitDataBlock"));

  if (!rddb){
    MSG("NtpSR",Msg::kWarning) 
      << "No RawDigitDataBlock Found, will return the default value"<<endl;
    return ratio;
  }

  int firstbucket = 999999999; //the first time bucket in the current event
  int firstbucket_corr = 999999999; //corrected firstbucket after removing early digits
 
  CandStripHandleItr stripItr(cndevt->GetDaughterIterator());
  while (CandStripHandle* strip=dynamic_cast<CandStripHandle*>(stripItr())){
    CandDigitHandleItr digitItr(strip->GetDaughterIterator());
    while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(digitItr())){
      const RawDigit* rd=dynamic_cast<const RawDigit*>(rddb->At(digit->GetRawDigitIndex()));
      if(rd->GetTDC()<firstbucket) firstbucket = rd->GetTDC();
      if(rd->GetTDC()<firstbucket_corr && rd->GetADC()>0.01*cndevt->GetCharge(CalStripType::kSigLin)){//get rid of early digit(s) with low ph (less than 1% of total ADC)
        firstbucket_corr = rd->GetTDC();
      }
    }
  }
  if (firstbucket_corr == 999999999) firstbucket_corr = firstbucket;

  float totadc = 0;  //total ADC in the current event
  float lateadc = 0; //total ADC after first 3 time buckets
  CandStripHandleItr stripItr2(cndevt->GetDaughterIterator());
  while (CandStripHandle* strip=dynamic_cast<CandStripHandle*>(stripItr2())){
    CandDigitHandleItr digitItr(strip->GetDaughterIterator());
    while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(digitItr())){
      const RawDigit* rd=dynamic_cast<const RawDigit*>(rddb->At(digit->GetRawDigitIndex()));
      totadc += rd->GetADC();
      if (rd->GetTDC()>=firstbucket_corr+3){
        lateadc += rd->GetADC();
      }      
    }
  }

  if (totadc) ratio = lateadc/totadc;

  return ratio;
}

//......................................................................
float NtpSRBleachFiller::GetstraightPHFraction(const CandEventHandle *evt,
                                               const CandRecord *crec){
  //
  // Purpose: Calculate the fraction of hits in an event that have previous
  //          activity in them.
  //
  // Arguments: CandEventHandle ptr, CandRecord ptr.
  // 
  // Return: straight ph fraction. -1 if error.
  //
  // Authors: Tom Osiecki (osiecki@mail.hep.utexas.edu)
  //
  // Creation date: 11/05 
  //

  if(!evt) {
    MSG("NtpSR",Msg::kWarning) << "No CandEventHandle Found in GetstraightPHFraction" << endl;
    return -1.0;
  }

  if(!crec) {
    MSG("NtpSR",Msg::kWarning) << "No CandRecord Found in GetstraightPHFraction" << endl;
    return -1.0;
  }

  const CandStripListHandle *striplisthandle = dynamic_cast <const CandStripListHandle*> 
    (crec->FindCandHandle("CandStripListHandle"));

  const CandStripListHandle *striplisthandle2 = dynamic_cast <const CandStripListHandle*> 
    (crec->FindCandHandle("CandStripListHandle"));

  if(!striplisthandle) return -1.0;

  TIter stripItr(striplisthandle->GetDaughterIterator());
  TIter stripItr2(striplisthandle2->GetDaughterIterator());
  
  int totalstrips = striplisthandle->GetNDaughters();      
  int *striparray = new int[totalstrips];
  std::map<Int_t, Int_t> uidmap;

  for(int zeroout = 0; zeroout < totalstrips; zeroout++) {
    striparray[zeroout] = -1;
  }
 
  Int_t index = 0;
  while(CandStripHandle* strip = dynamic_cast<CandStripHandle*> (stripItr())) {
    Int_t uid = strip -> GetUidInt();
    
    stripItr2.Reset();
    while(CandStripHandle* strip2 = dynamic_cast<CandStripHandle*> (stripItr2())) {
      Int_t uid2 = strip2 -> GetUidInt();
          
      if(uid != uid2 && 
         strip -> GetStrip() == strip2 -> GetStrip() && 
         strip -> GetPlane() == strip2 -> GetPlane()) {
      
        if((strip -> GetTime() - strip2 -> GetTime())*1.0e9 > 0.0) {
          striparray[index]++; 
          uidmap[uid] = index;
        }
        
      } // Don't compare the same strip to itself.
    }
    index++;
  }
  
  Int_t nstp = 0, nstpevt = 0;
  CandStripHandleItr evtstripItr(evt->GetDaughterIterator());
  while (CandStripHandle* strip = dynamic_cast<CandStripHandle*> (evtstripItr())) {
    nstpevt++;
    if(striparray[uidmap[strip -> GetUidInt()]] > -1.0) nstp++;
  }

  Float_t ratio = (float) nstp / nstpevt;
  MSG("NtpSR",Msg::kDebug) << "GetstraightPHFraction ratio = " << ratio << endl;

  delete [] striparray; striparray = 0;

  return ratio;
}

---