CalHit.cxx

Go to the documentation of this file.

#include <sstream>
#include "RawData/RawDigit.h"
#include "CalHit/CalHit.h"
#include "JobControl/Exception.h"
#include "Geometry/Geometry.h"
#include "Geometry/PlaneGeo.h"
#include "Geometry/CellGeo.h"

using namespace calhit;

MODULE_DECL(CalHit);

CalHit::CalHit(const char* version) : 
jobc::Module("CalHit"),
fInputFolder("."),
fOutputFolder("Hits"),
fDoCalib(1),
fWarnedPE(false),
fWarnedPEAC(false),
fWarnedMIP(false),
fWarnedNS(false),
fWarnedTT0(false),
fWarnedTTrans(false)
{

   this->SetCfgVersion(version);
}

CalHit::~CalHit()
{

}

jobc::Result CalHit::Reco(edm::EventHandle& evt)
{

  // base class Reco method.  

   if (jobc::kFailed==GetData(evt)) {
           std::cerr <<"CalHit::GetData failed.  Returning"<<std::endl;
           return jobc::kFailed;
   }

   if (jobc::kFailed==AddToData(evt)) {
           std::cerr <<"CalHit::AddToData failed.  Returning"<<std::endl;
           return jobc::kFailed;
   }

   if (jobc::kFailed==SaveData(evt)) {
           std::cerr <<"CalHit::SaveData failed.  Returning"<<std::endl;
           return jobc::kFailed;
   }

 
   return jobc::kPassed;

}


jobc::Result CalHit::GetData(edm::EventHandle& evt)
{
      // clear the local copy
      fCellVec.clear();
      fCellInput.clear();

       // get the data from the appropriate place
      if (DoCalib(recobase::kRAW)) {

         // create cal hits, and copy over raw data
             if (Raw(evt)==jobc::kFailed) return jobc::kFailed;

      } else {

         // start with existing cal hits
             if (Read(evt)==jobc::kFailed) return jobc::kFailed;

      }

      return jobc::kPassed;
}

jobc::Result CalHit::Raw(edm::EventHandle& evt)
{
    // get raw data
    
    std::vector<const rawdata::RawDigit*> diglist;
    try {
       evt.Raw().Get(fInputFolder.c_str(),diglist);
    } catch (edm::Exception e) {
       std::cerr << "Error retrieving digit list, while looking for raw "<<
         "digits in CalHit::Raw()"<<std::endl;
       return jobc::kFailed;
    }
  
    // loop over all raw digits, and create basic calibrated hits in local
    // record
    std::vector<const rawdata::RawDigit*>::iterator iter=diglist.begin();
    // get geometry N.S.
    short fDetGeom=1;
    fGeo = &geo::Geometry::Instance(evt.Header().Run(), fDetGeom);
    
    while (iter!=diglist.end())
    {
       const rawdata::RawDigit *rd=*iter;
       recobase::CellHit* ch=new recobase::CellHit(rd->Channel(),
                                                      rd->ADC(0),rd->TDC(0));
        ch->SetCell (rd->GetCell());
        ch->SetPlane(rd->GetPlane());
        
        // N.S 08/25/08
        double posit[3];

        fGeo->Plane(rd->GetPlane()).Cell(rd->GetCell()).GetCenter(posit,0.0);
        if(fGeo->Plane(rd->GetPlane()).View()==geo::kY) {
        ch->SetTCpos(posit[1]);
        ch->SetPLVC(2);
        }
        else{
        ch->SetTCpos(posit[0]);
        ch->SetPLVC(1);
        }
        ch->SetZCPos(posit[2]);
        // End N.S. 08/25/08
                
        fCellVec.push_back(ch);

       // if we'll be doing an attenuation correction, collect
       // necessary information here
       if (DoCalib(recobase::kPEAC)) AddAttCorr(ch);
       iter++;
    }
    
    return jobc::kPassed;

}
  
jobc::Result CalHit::Read(edm::EventHandle& evt)
{
    // get existing CalHit's
    
    try {
       evt.Reco().Get(fInputFolder.c_str(),fCellInput);
    } catch (edm::Exception e) {
       std::cerr << "Error retrieving CellHit list in "<<fInputFolder.c_str()
               <<std::endl;
       return jobc::kFailed;
    }
  

  // if we will need attenuation correction, accumulate necessary info
  // here
    std::vector<const recobase::CellHit*>::iterator iter=fCellInput.begin();
    std::cerr<<"Found "<<fCellInput.size()<<" cell hits in input"<<std::endl;
    while (iter!=fCellInput.end())
    {
       const recobase::CellHit* oldch=*iter;
       recobase::CellHit* ch=new recobase::CellHit(*oldch);
       fCellVec.push_back(ch);
       if (DoCalib(recobase::kPEAC)) 
       {
          AddAttCorr(ch);
       }
       iter++;
    }

    return jobc::kPassed;
}

jobc::Result CalHit::AddToData(edm::EventHandle& evt)
{

   std::vector<recobase::CellHit*>::iterator iter=fCellVec.begin();
   while (iter!=fCellVec.end())
   {
       recobase::CellHit* ch=*iter;

       // create Photo Electron equivalent ADC
       if (DoCalib(recobase::kPE))
                if (ADC2PEHit(ch)==jobc::kFailed) return jobc::kFailed;

       //... with Attenuation correction
       if (DoCalib(recobase::kPEAC))
                if (ADC2PEACHit(ch)==jobc::kFailed) return jobc::kFailed;

       //... in MIP units
       if (DoCalib(recobase::kMIP))
                if (ADC2MIPHit(ch)==jobc::kFailed) return jobc::kFailed;

       // tdc in ns
       if (DoCalib(recobase::kTNS))
                if (TDC2NSHit(ch)==jobc::kFailed) return jobc::kFailed;

       // tdc in ns, with transit time correction
       if (DoCalib(recobase::kTTrans))
                if (TDC2TTransHit(ch)==jobc::kFailed) return jobc::kFailed;

       // tdc in ns, with transit time correction, from a T0
       if (DoCalib(recobase::kTT0))
                if (TDC2TT0Hit(ch)==jobc::kFailed) return jobc::kFailed;

       iter++;
   }

   return jobc::kPassed;

}

jobc::Result CalHit::SaveData(edm::EventHandle& evt)
{
   if (evt.Reco().GetFolder(fOutputFolder.c_str())) {
       std::cerr << "Output Folder reco/"<<fOutputFolder.c_str()<<
         " already exists"<<std::endl;
       return jobc::kFailed;
   }

   TFolder* folder=evt.Reco().MakeFolder(fOutputFolder.c_str());

   unsigned int nsave=evt.Reco().PutVector(fCellVec,fOutputFolder.c_str());
   if (nsave!=fCellVec.size())
   {
       std::cerr <<"CalHit::SaveData: wrong number CellHits saved: "<<
          nsave<<" instead of expected "<<fCellVec.size()<<std::endl;
       return jobc::kFailed;
   }
   float pe,peac;
   fCellVec[0]->PE(pe);
   fCellVec[0]->PEAC(peac);
   std::cerr << "saved "<<nsave<<" CellHits.  first pe="<<pe<<", atcor="<<peac<<
         std::endl;

   return jobc::kPassed;
       
}

bool CalHit::DoCalib(recobase::CellHit_t level, int set)
{
    // if set<0, then caller is asking what calibrations are done
    //    set=0, caller wants to turn off a calibration
    //    set>0, caller wants to turn on a calibration

     if (set < 0 ) {
         return (fDoCalib & (1<<level) );
     } else if (set == 0 ) {
         fDoCalib = fDoCalib&=(~(1<<level));
         return false;
     } else  {
         fDoCalib = fDoCalib|=(1<<level);
         return true;
     }
     
}


void CalHit::AddAttCorr(recobase::CellHit* ch)
{
   // add this cell to attenuation correction

   // but a dummy for the time being...


}

jobc::Result CalHit::ADC2PEHit(recobase::CellHit* ch)
{
    // Boneheaded base clase
    if (!fWarnedPE) {
       std::cerr<<"\nWarning!  Using baseclass CalHit::ADC2PEHit. "<<
                   " Do not use this for Physics!"<<std::endl;
       fWarnedPE=true; 
    }  


    unsigned short adc;
    if (!ch->ADC(adc)) {
       std::cerr<<"CalHit::ADC2PEHit(): raw adc not set in cell hit for channel "<<
                  ch->Channel()<<std::endl;
       return jobc::kFailed;
    }
    ch->SetPE(adc/4.095);// dummy value

    return jobc::kPassed;

} // ADC2PEHit

jobc::Result CalHit::ADC2PEACHit(recobase::CellHit* ch)
{
    // Boneheaded base clase - represents attenuation-corrected PEs

    if (!fWarnedPEAC) {
       std::cerr<<"\nWarning!  Using baseclass CalHit::ADC2PEACHit. "<<
                   " Do not use this for Physics!"<<std::endl;
       fWarnedPEAC=true; 
    }  


  // require PE number before applying attenuation correction
    float pe;
    if (!ch->PE(pe)) {
           if (!ADC2PEHit(ch)) return jobc::kFailed;
           ch->PE(pe);
    }
    ch->SetPEAC(pe/0.55);

    return jobc::kPassed;

} // ADC2PEACHit

jobc::Result CalHit::ADC2MIPHit(recobase::CellHit* ch)
{
    // Boneheaded base clase - represents MIP-corrected ph

    if (!fWarnedMIP) {
       std::cerr<<"\nWarning!  Using baseclass CalHit::ADC2MIPHit. "<<
                   " Do not use this for Physics!"<<std::endl;
       fWarnedMIP=true; 
    }  


  // require PE attcor'd number before converting to MIPs
    float pe;
    if (!ch->PEAC(pe)) {
           if (!ADC2PEACHit(ch)) return jobc::kFailed;
           ch->PEAC(pe);
    }
    ch->SetMIP(pe/30.);

    return jobc::kPassed;

} // ADC2MIPHit

jobc::Result CalHit::TDC2NSHit(recobase::CellHit* ch)
{
    // Boneheaded base clase - represents time in ns
    if (!fWarnedNS) {
       std::cerr<<"\nWarning!  Using baseclass CalHit::TDC2NSHit. "<<
                   " Do not use this for Physics!"<<std::endl;
       fWarnedNS=true; 
    }  


    unsigned short tdc;
    if (!ch->TDC(tdc)) {
       std::cerr<<"CalHit::TDC2NSHit(): raw tdc not set in cell hit  for channel "<<
                  ch->Channel()<<std::endl;
       return jobc::kFailed;
    }
    ch->SetTNS(tdc*1.5e3); 

    return jobc::kPassed;

} // TDC2NSHit

jobc::Result CalHit::TDC2TT0Hit(recobase::CellHit* ch)
{
    // Boneheaded base clase - represents hit times corrected for global event
    // t0

    if (!fWarnedTT0) {
       std::cerr<<"\nWarning!  Using baseclass CalHit::TDC2TT0Hit. "<<
                   " Do not use this for Physics!"<<std::endl;
       fWarnedTT0=true; 
    }  


  // require NS number before applying t0
    double ns;
    if (!ch->TNS(ns)) {
           if (!TDC2NSHit(ch)) return jobc::kFailed;
           ch->TNS(ns); // should be there now...
    }
    ch->SetTT0(ns+24.6);

    return jobc::kPassed;

} // TDC2TT0Hit

jobc::Result CalHit::TDC2TTransHit(recobase::CellHit* ch)
{
    // Boneheaded base clase - represents transit-time corrected hits

    if (!fWarnedTTrans) {
       std::cerr<<"\nWarning!  Using baseclass CalHit::TDC2TTransHit. "<<
                   " Do not use this for Physics!"<<std::endl;
       fWarnedTTrans=true; 
    }  

  // require t0-corrected time before handling transit time
    double tt0;
    if (!ch->TT0(tt0)) {
           if (!TDC2TT0Hit(ch)) return jobc::kFailed; // try setting ns
           ch->TT0(tt0); // should be there now...
    }
    ch->SetTTrans(tt0-7.5*5); 

    return jobc::kPassed;

} // TDC2TTransHit

void CalHit::Update(const cfg::Config& c)
{
  c("Input"). Get(fInputFolder);
  c("Output"). Get(fOutputFolder);



    if (fInputFolder==fOutputFolder) {
       std::cerr<<"CalHit::CalHit() - output folder ("<<fOutputFolder.c_str()<<
          ") is not allowed to equal input folder ("<<fInputFolder.c_str()
                 <<")"<< std::endl;
       assert(0);
    }

   std::cerr<<"Set Input Folder ="<<fInputFolder.c_str()<<std::endl;
   std::cerr<<"Set Output Folder ="<<fOutputFolder.c_str()<<std::endl;


  // zero out calibrations
  fDoCalib=0;

  std::string calibs;
  c("Calibs"). Get(calibs);

  // set calibrations according to the string
  std::istringstream icalibs(calibs);
  while (icalibs.good())
  {
     std::string tmp;
     icalibs>>tmp;
     if (tmp=="RAW") DoCalib(recobase::kRAW,1);
     else if (tmp=="PE") DoCalib(recobase::kPE,1);
     else if (tmp=="PEAC")  DoCalib(recobase::kPEAC,1);
     else if (tmp=="MIP")  DoCalib(recobase::kMIP,1);
     else if (tmp=="TNS")  DoCalib(recobase::kTNS,1);
     else if (tmp=="TT0")  DoCalib(recobase::kTT0,1);
     else if (tmp=="TTRANS")  DoCalib(recobase::kTTrans,1);
     else {
         std::cerr<<"Unknown calibration requested: "<<tmp.c_str()<<std::endl;
         abort();
     }


  } // parse calibration string

  PrintCalibs();
  

}
void CalHit::PrintCalibs()
{
   std::cerr<<"Calibrations Requested:"<<std::endl;

   std::cerr<<" Starting from Raw Data: " ;
     if (DoCalib(recobase::kRAW)) std::cerr<<" Yes"<<std::endl; else
                         std::cerr<<" NO"<<std::endl;

   std::cerr<<" PE Units: ";
     if (DoCalib(recobase::kPE)) std::cerr<<" Yes"<<std::endl; else
                         std::cerr<<" NO"<<std::endl;

   std::cerr<<" Attenuation Correction (to PEs): ";
     if (DoCalib(recobase::kPEAC)) std::cerr<<" Yes"<<std::endl; else
                         std::cerr<<" NO"<<std::endl;

   std::cerr<<" MIP Units: ";
     if (DoCalib(recobase::kMIP)) std::cerr<<" Yes"<<std::endl; else
                         std::cerr<<" NO"<<std::endl;

   std::cerr<<" Time in ns units: " ;
     if (DoCalib(recobase::kTNS)) std::cerr<<" Yes"<<std::endl; else
                         std::cerr<<" NO"<<std::endl;

   std::cerr<<" T-zero corrected time in ns: " ;
     if (DoCalib(recobase::kTT0)) std::cerr<<" Yes"<<std::endl; else
                         std::cerr<<" NO"<<std::endl;

   std::cerr<<" Transit time corrected time in ns: " ;
     if (DoCalib(recobase::kTTrans)) std::cerr<<" Yes"<<std::endl; else
                         std::cerr<<" NO"<<std::endl;
}

//} // namespace calhit

---