IceXtalk.cxx

Go to the documentation of this file.

 
#include "IceXtalk.h"
#include "Riostream.h"
 
ClassImp(IceXtalk); // Class implementation to enable ROOT I/O

IceXtalk::IceXtalk(const char* name,const char* title) : TTask(name,title)
{
 
 fCalfile=0;
 fOmdb=0;
 fPmin=0.5;
 fPe=1;
}

IceXtalk::~IceXtalk()
{
 
 if (fCalfile)
 {
  delete fCalfile;
  fCalfile=0;
 }
}

void IceXtalk::SetOMdbase(NcObjMatrix* omdb)
{
 
 fOmdb=omdb;
}

void IceXtalk::SetCalibFile(TString name)
{
 
 // Expand the path name of the provided input file
 name=gSystem->ExpandPathName(name.Data());
 
 if (fCalfile)
 {
  delete fCalfile;
  fCalfile=0;
 }
 fCalfile=new TFile(name.Data());
 fOmdb=0;
}

void IceXtalk::SetMinProb(Float_t pmin)
{
 
 fPmin=pmin;
}

void IceXtalk::SetXtalkPE(Float_t pe)
{
 
 fPe=pe;
}

void IceXtalk::Exec(Option_t* opt)
{
 
 TString name=opt;
 NcJob* parent=(NcJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
 
 if (!parent) return;
 
 IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent");
 if (!evt) return;
 
 // Only process accepted events
 NcDevice* seldev=(NcDevice*)evt->GetDevice("NcEventSelector");
 if (seldev)
 {
  if (seldev->GetSignal("Select") < 0.1) return;
 }
 
 NcSignal* daq=(NcSignal*)evt->GetDevice("Daq");
 Int_t mudaq=0;
 if (daq) mudaq=int(daq->GetSignal("Muon"));
 
 // This cross talk correction processor is only for MuDaq data 
 if (!mudaq) return;
 
 if (!fOmdb && fCalfile)
 {
  fOmdb=(NcObjMatrix*)fCalfile->Get("MuDaq-OMDBASE");
  // Next statement for compatibility with old calibration file format
  if (!fOmdb) fOmdb=(NcObjMatrix*)fCalfile->Get("Cal-OMDBASE");
 }
 
 if (!fOmdb) return;
 
 // Storage of the used parameters in the IceXtalk device
 NcDevice params;
 params.SetNameTitle("IceXtalk","IceXtalk processor parameters");
 params.SetSlotName("Pmin",1);
 params.SetSlotName("Pe",2);
 
 params.SetSignal(fPmin,1);
 params.SetSignal(fPe,2);
 
 evt->AddDevice(params);
 
 // All Amanda OMs with a signal
 TObjArray* mods=evt->GetDevices("IceAOM");
 if (!mods) return;
 Int_t nmods=mods->GetEntries();
 if (!nmods) return;
 
 TObjArray xhits; // Array with pointers to Xtalk hits to be corrected
 
 IceAOM* omt=0; // Transmitter OM
 IceAOM* omr=0; // Receiver OM
 TF1* fxtalk=0; // The Xtalk probability function
 Int_t idtrans=0;
 Int_t idrec=0;
 Float_t adct=0,let=0;
 Float_t p=0;
 Float_t adcr=0,ler=0;
 Int_t ibad=0;
 Float_t cpar=0,bpar=0,dlemin=0,dlemax=0,test=0;
 Float_t dle=0;
 Float_t sigcor=0;
 for (Int_t imod=0; imod<nmods; imod++)
 {
  omt=(IceAOM*)mods->At(imod);
  if (!omt) continue;
  idtrans=omt->GetUniqueID();
 
  // Also take bad transmitter modules into account
  ibad=omt->GetDeadValue("ADC");
  if (ibad) omt->SetAlive("ADC");
 
  // Check for corresponding receiver modules  
  for (Int_t jmod=0; jmod<nmods; jmod++)
  {
   // No Xtalk from a module to itself 
   if (jmod==imod) continue;
 
   omr=(IceAOM*)mods->At(jmod);
   if (!omr) continue;
   idrec=omr->GetUniqueID();
 
   fxtalk=(TF1*)fOmdb->GetObject(idtrans,idrec+1);
   if (!fxtalk) continue;
 
   // Determine Xtalk probability for each transmitter hit
   for (Int_t ithit=1; ithit<=omt->GetNhits(); ithit++)
   {
    NcSignal* st=omt->GetHit(ithit);
    if (!st) continue;
 
    // Eliminate possible previous Xtalk correction
    sigcor=st->GetOffset("ADC");
    st->AddSignal(sigcor,"ADC");
    st->ResetOffset("ADC");
    adct=st->GetSignal("ADC",-4); // Get uncalibrated amplitude
    dlemin=fxtalk->GetParameter("dLE-min");
    dlemax=fxtalk->GetParameter("dLE-max");
    // Protection against overflow in Xtalk prob. function
    cpar=fxtalk->GetParameter("C");
    bpar=fxtalk->GetParameter("B");
    test=(adct-cpar)/bpar;
    if (test>100)
    {
     p=1;
    }
    else if (test<-100)
    {
     p=0;
    }
    else
    {
     p=fxtalk->Eval(adct);
    }
    if (p<fPmin) continue;
 
    // Xtalk flagging for each receiver hit
    for (Int_t irhit=1; irhit<=omr->GetNhits(); irhit++)
    {
     NcSignal* sr=omr->GetHit(irhit);
     if (!sr) continue;
 
     // Check calibrated LE time differences
     if (!ibad)
     {
      let=st->GetSignal("LE",4);
      ler=sr->GetSignal("LE",4);
      dle=ler-let;
      if (dle<dlemin || dle>dlemax) continue;
     }
     // Register this receiver hit to be corrected
     xhits.Add(sr); 
    }
   } // end of transmitter hit loop
  } // end of receiver OM loop
 
  // Restore the original transmitter dead flag value
  if (ibad) omt->SetDead("ADC");
 } // end of transmitter OM loop
 
 // Perform the Xtalk signal correction to the registered receiver hits
 // The correction value is stored in the signal data as an offset value
 for (Int_t ix=0; ix<xhits.GetEntries(); ix++)
 {
  NcSignal* sx=(NcSignal*)xhits.At(ix);
  if (!sx) continue;
 
  // Eliminate possible previous Xtalk correction
  sigcor=sx->GetOffset("ADC");
  adcr=sx->GetSignal("ADC")+sigcor;
  sigcor=fPe;
  // If stored ADC data is un-calibrated, convert fPe to raw ADC value 
  NcSignal* parent=(NcSignal*)sx->GetDevice();
  if (parent)
  {
   if (parent->GetCalFunction("ADC"));
   {
    idrec=parent->GetUniqueID();
    omr=(IceAOM*)fOmdb->GetObject(idrec,1);
    TF1* fdecal=0;
    if (omr) fdecal=omr->GetDecalFunction("ADC");
    if (fdecal) sigcor=fdecal->Eval(fPe);
   }
  }
  adcr=adcr-sigcor;
  sx->SetSignal(adcr,"ADC");
  sx->SetOffset(sigcor,"ADC");
 }
}