IceCleanHits.cxx

Go to the documentation of this file.

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

IceCleanHits::IceCleanHits(const char* name,const char* title) : TTask(name,title)
{
 
 fEvt=0;
 fAdcminAM=0.3;
 fAdcmaxAM=999999;
 fAdcminAT=0.3;
 fAdcmaxAT=999999;
 fAdcminT=0;
 fAdcmaxT=999999;
 fAdcminIC=0;
 fAdcmaxIC=999999;
 fAdcminDC=0;
 fAdcmaxDC=999999;
 fTotminAEM=125;
 fTotmaxAEM=2000;
 fTotminAOM=20;
 fTotmaxAOM=2000;
 fTotminAET=125;
 fTotmaxAET=2000;
 fTotminAOT=20;
 fTotmaxAOT=2000;
 fTotminT=0;
 fTotmaxT=999999;
 fTotminIC=0;
 fTotmaxIC=999999;
 fTotminDC=0;
 fTotmaxDC=999999;
 fRmaxA=70;
 fDtmaxA=500;
 fRmaxT=150;
 fDtmaxT=1000;
 fRmaxIC=150;
 fDtmaxIC=1000;
 fRmaxDC=75;
 fDtmaxDC=500;
 fTnamAM="main";
 fTtimAM=-1;
 fTnamAT="main";
 fTtimAT=-1;
 fTnamT="user";
 fSlcT=0;
 fMultT=1;
 fTtimT=-1;
 fTnamIC="user";
 fSlcIC=0;
 fMultIC=1;
 fTtimIC=-1;
 fTnamDC="user";
 fSlcDC=0;
 fMultDC=1;
 fTtimDC=-1;
 fTwinAM=2250;
 fTwinAT=3000;
 fTwinT=999999;
 fTwinIC=3000;
 fTwinDC=1000;
}

IceCleanHits::~IceCleanHits()
{
}

void IceCleanHits::SetAdcRange(Float_t min,Float_t max,TString s)
{
 
 if (s=="MuDaq")
 {
  fAdcminAM=min;
  fAdcmaxAM=max;
 }
 if (s=="TWRDaq")
 {
  fAdcminAT=min;
  fAdcmaxAT=max;
 }
 if (s=="TDOM")
 {
  fAdcminT=min;
  fAdcmaxT=max;
 }
 if (s=="ICDOM")
 {
  fAdcminIC=min;
  fAdcmaxIC=max;
 }
 if (s=="DCDOM")
 {
  fAdcminDC=min;
  fAdcmaxDC=max;
 }
}

void IceCleanHits::SetTotRange(Float_t min,Float_t max,TString s)
{
 
 if (s=="AEM")
 {
  fTotminAEM=min;
  fTotmaxAEM=max;
 }
 if (s=="AOM")
 {
  fTotminAOM=min;
  fTotmaxAOM=max;
 }
 if (s=="AET")
 {
  fTotminAET=min;
  fTotmaxAET=max;
 }
 if (s=="AOT")
 {
  fTotminAOT=min;
  fTotmaxAOT=max;
 }
 if (s=="TDOM")
 {
  fTotminT=min;
  fTotmaxT=max;
 }
 if (s=="ICDOM")
 {
  fTotminIC=min;
  fTotmaxIC=max;
 }
 if (s=="DCDOM")
 {
  fTotminDC=min;
  fTotmaxDC=max;
 }
}

void IceCleanHits::SetIsolation(Float_t rmax,Float_t dtmax,TString s)
{
 
 if (s=="AOM")
 {
  if (rmax >= 0) fRmaxA=rmax;
  if (dtmax >= 0) fDtmaxA=dtmax;
 }
 if (s=="TDOM")
 {
  if (rmax >= 0) fRmaxT=rmax;
  if (dtmax >= 0) fDtmaxT=dtmax;
 }
 if (s=="ICDOM")
 {
  if (rmax >= 0) fRmaxIC=rmax;
  if (dtmax >= 0) fDtmaxIC=dtmax;
 }
 if (s=="DCDOM")
 {
  if (rmax >= 0) fRmaxDC=rmax;
  if (dtmax >= 0) fDtmaxDC=dtmax;
 }
}

void IceCleanHits::SetTwindow(Float_t dtmax,TString s)
{
 
 if (s=="MuDaq") fTwinAM=dtmax;
 if (s=="TWRDaq") fTwinAT=dtmax;
 if (s=="TDOM") fTwinT=dtmax;
 if (s=="ICDOM") fTwinIC=dtmax;
 if (s=="DCDOM") fTwinDC=dtmax;
}

void IceCleanHits::SetTname(TString name,TString s)
{
 
 if (s=="MuDaq")
 {
  fTnamAM=name;
  fTtimAM=0;
  if (name=="Average") fTtimAM=-1;
 }
 if (s=="TWRDaq")
 {
  fTnamAT=name;
  fTtimAT=0;
  if (name=="Average") fTtimAT=-1;
 }
 if (s=="TDOM")
 {
  fTnamT=name;
  fTtimT=0;
  if (name=="Average") fTtimT=-1;
 }
 if (s=="ICDOM")
 {
  fTnamIC=name;
  fTtimIC=0;
  if (name=="Average") fTtimIC=-1;
 }
 if (s=="DCDOM")
 {
  fTnamDC=name;
  fTtimDC=0;
  if (name=="Average") fTtimDC=-1;
 }
}

void IceCleanHits::SetTtime(Float_t t,TString s,Int_t slc,Int_t mult)
{
 
 if (s=="MuDaq")
 {
  fTtimAM=t;
  fTnamAM="user";
 }
 if (s=="TWRDaq")
 {
  fTtimAT=t;
  fTnamAT="user";
 }
 if (s=="TDOM")
 {
  fTtimT=t;
  fSlcT=slc;
  fMultT=mult;
  fTnamT="user";
 }
 if (s=="ICDOM")
 {
  fTtimIC=t;
  fSlcIC=slc;
  fMultIC=mult;
  fTnamIC="user";
 }
 if (s=="DCDOM")
 {
  fTtimDC=t;
  fSlcDC=slc;
  fMultDC=mult;
  fTnamDC="user";
 }
}

void IceCleanHits::Exec(Option_t* opt)
{
 
 TString name=opt;
 NcJob* parent=(NcJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
 
 if (!parent) return;
 
 fEvt=(IceEvent*)parent->GetObject("IceEvent");
 if (!fEvt) return;
 
 // Only process accepted events
 NcDevice* seldev=(NcDevice*)fEvt->GetDevice("NcEventSelector");
 if (seldev)
 {
  if (seldev->GetSignal("Select") < 0.1) return;
 }
 
 // Storage of the used parameters in the IceCleanHits device
 NcDevice params;
 params.SetNameTitle("IceCleanHits","IceCleanHits processor parameters");
 params.AddNamedSlot("AdcminAM");
 params.AddNamedSlot("AdcmaxAM");
 params.AddNamedSlot("AdcminAT");
 params.AddNamedSlot("AdcmaxAT");
 params.AddNamedSlot("AdcminT");
 params.AddNamedSlot("AdcmaxT");
 params.AddNamedSlot("AdcminIC");
 params.AddNamedSlot("AdcmaxIC");
 params.AddNamedSlot("AdcminDC");
 params.AddNamedSlot("AdcmaxDC");
 params.AddNamedSlot("TotminAEM");
 params.AddNamedSlot("TotmaxAEM");
 params.AddNamedSlot("TotminAOM");
 params.AddNamedSlot("TotmaxAOM");
 params.AddNamedSlot("TotminAET");
 params.AddNamedSlot("TotmaxAET");
 params.AddNamedSlot("TotminAOT");
 params.AddNamedSlot("TotmaxAOT");
 params.AddNamedSlot("TotminT");
 params.AddNamedSlot("TotmaxT");
 params.AddNamedSlot("TotminIC");
 params.AddNamedSlot("TotmaxIC");
 params.AddNamedSlot("TotminDC");
 params.AddNamedSlot("TotmaxDC");
 params.AddNamedSlot("RmaxA");
 params.AddNamedSlot("DtmaxA");
 params.AddNamedSlot("RmaxT");
 params.AddNamedSlot("DtmaxT");
 params.AddNamedSlot("RmaxIC");
 params.AddNamedSlot("DtmaxIC");
 params.AddNamedSlot("RmaxDC");
 params.AddNamedSlot("DtmaxDC");
 params.AddNamedSlot("TwinAM");
 params.AddNamedSlot("TtimAM");
 params.AddNamedSlot("TwinAT");
 params.AddNamedSlot("TtimAT");
 params.AddNamedSlot("TwinT");
 params.AddNamedSlot("SlcT");
 params.AddNamedSlot("MultT");
 params.AddNamedSlot("TtimT");
 params.AddNamedSlot("TwinIC");
 params.AddNamedSlot("SlcIC");
 params.AddNamedSlot("MultIC");
 params.AddNamedSlot("TtimIC");
 params.AddNamedSlot("TwinDC");
 params.AddNamedSlot("SlcDC");
 params.AddNamedSlot("MultDC");
 params.AddNamedSlot("TtimDC");
 
 params.SetSignal(fAdcminAM,"AdcminAM");
 params.SetSignal(fAdcmaxAM,"AdcmaxAM");
 params.SetSignal(fAdcminAT,"AdcminAT");
 params.SetSignal(fAdcmaxAT,"AdcmaxAT");
 params.SetSignal(fAdcminT,"AdcminT");
 params.SetSignal(fAdcmaxT,"AdcmaxT");
 params.SetSignal(fAdcminIC,"AdcminIC");
 params.SetSignal(fAdcmaxIC,"AdcmaxIC");
 params.SetSignal(fAdcminDC,"AdcminDC");
 params.SetSignal(fAdcmaxDC,"AdcmaxDC");
 params.SetSignal(fTotminAEM,"TotminAEM");
 params.SetSignal(fTotmaxAEM,"TotmaxAEM");
 params.SetSignal(fTotminAOM,"TotminAOM");
 params.SetSignal(fTotmaxAOM,"TotmaxAOM");
 params.SetSignal(fTotminAET,"TotminAET");
 params.SetSignal(fTotmaxAET,"TotmaxAET");
 params.SetSignal(fTotminAOT,"TotminAOT");
 params.SetSignal(fTotmaxAOT,"TotmaxAOT");
 params.SetSignal(fTotminT,"TotminT");
 params.SetSignal(fTotmaxT,"TotmaxT");
 params.SetSignal(fTotminIC,"TotminIC");
 params.SetSignal(fTotmaxIC,"TotmaxIC");
 params.SetSignal(fTotminDC,"TotminDC");
 params.SetSignal(fTotmaxDC,"TotmaxDC");
 params.SetSignal(fRmaxA,"RmaxA");
 params.SetSignal(fDtmaxA,"DtmaxA");
 params.SetSignal(fRmaxT,"RmaxT");
 params.SetSignal(fDtmaxT,"DtmaxT");
 params.SetSignal(fRmaxIC,"RmaxIC");
 params.SetSignal(fDtmaxIC,"DtmaxIC");
 params.SetSignal(fRmaxDC,"RmaxDC");
 params.SetSignal(fDtmaxDC,"DtmaxDC");
 params.SetSignal(fTwinAM,"TwinAM");
 params.SetSignal(fTtimAM,"TtimAM");
 params.SetSignal(fTwinAT,"TwinAT");
 params.SetSignal(fTtimAT,"TtimAT");
 params.SetSignal(fTwinT,"TwinT");
 params.SetSignal(fSlcT,"SlcT");
 params.SetSignal(fMultT,"MultT");
 params.SetSignal(fTtimT,"TtimT");
 params.SetSignal(fTwinIC,"TwinIC");
 params.SetSignal(fSlcIC,"SlcIC");
 params.SetSignal(fMultIC,"MultIC");
 params.SetSignal(fTtimIC,"TtimIC");
 params.SetSignal(fTwinDC,"TwinDC");
 params.SetSignal(fSlcDC,"SlcDC");
 params.SetSignal(fMultDC,"MultDC");
 params.SetSignal(fTtimDC,"TtimDC");
 
 fEvt->AddDevice(params);
 
 Amanda();
 IceTop();
 ICdoms();
 DCdoms();
}

void IceCleanHits::Amanda()
{
 
 NcDevice* daq=(NcDevice*)fEvt->GetDevice("Daq");
 if (!daq) return;
 
 if (daq->GetSignal("Muon")) MuDaq();
 if (daq->GetSignal("TWR")) TWRDaq();
}

void IceCleanHits::MuDaq()
{
 
 // All Amanda OMs with a signal
 TObjArray* aoms=fEvt->GetDevices("IceAOM");
 if (!aoms) return;
 
 // Local OM array with bad/dead OMs (as indicated via IceCalibrate) discarded
 TObjArray oms;
 IceAOM* omx=0;
 for (Int_t i=0; i<aoms->GetEntries(); i++)
 {
  omx=(IceAOM*)aoms->At(i);
  if (!omx) continue;
  if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
  oms.Add(omx);
 }
 
 // Determination of the selected trigger time
 if (fTnamAM != "user") // Trigger time derived from trigger data
 {
  fTtimAM=fEvt->GetTriggerTime(fTnamAM);
 }
 else // Trigger time from median of hit times
 {
  if (fTtimAM<0)
  {
   fTtimAM=fEvt->GetTriggerTime(fTnamAM,&oms);
 
   // Make the value negative to stay with median of hit times selection
   // for all following invokations.
   if (fTtimAM>0) fTtimAM*=-1.;
   if (!fTtimAM) fTtimAM=-0.001;
  }
 }
 
 // Store the trigger time in the parameters device
 NcDevice* params=(NcDevice*)fEvt->GetDevice("IceCleanHits");
 if (params) params->SetSignal(fTtimAM,"TtimAM");
 
 // Local array with clean hits
 TObjArray hits;
 Int_t omid=0;
 Int_t clean=1;
 NcSignal* sx=0;
 Float_t adc,le,tot;
 Int_t readout;
 for (Int_t iom=0; iom<oms.GetEntries(); iom++)
 {
  omx=(IceAOM*)oms.At(iom);
  if (!omx) continue;
  omid=omx->GetUniqueID();
  readout=int(omx->GetSignal("READOUT"));
  // General readout setting in case info was missing 
  if (!readout)
  {
   readout=1;
   if (omid>=303) readout=2; // Optical OMs
  }
  for (Int_t ih=1; ih<=omx->GetNhits(); ih++)
  {
   sx=omx->GetHit(ih);
   if (!sx) continue;
   adc=sx->GetSignal("ADC",7);
   le=sx->GetSignal("LE",7);
   tot=sx->GetSignal("TOT",7);
 
   // Remove hits with an ADC value outside the range
   if (adc<fAdcminAM || adc>fAdcmaxAM)
   {
    sx->SetDead("ADC");
    clean=0;
   }
   // Remove hits with a TOT value outside the range
   // Note : Different ranges for electrical and optical modules
   if (readout==1) // Electrical OMs
   {
    if (tot<fTotminAEM || tot>fTotmaxAEM)
    {
     sx->SetDead("TOT");
     clean=0;
    }
   }
   else // Optical OMs
   {
    if (tot<fTotminAOM || tot>fTotmaxAOM)
    {
     sx->SetDead("TOT");
     clean=0;
    }
   }
   // Remove hits that are outside the trigger time window.
   // Since the trigger time was determined from uncalibrated LE's
   // (to include cable length effects) the uncalibrated LE of each
   // hit should be used here as well. 
   le=sx->GetSignal("LE",-7);
   if (fabs(le-fTtimAM)>fTwinAM)
   {
     sx->SetDead("LE");
     clean=0;
   }
   // Store only the current clean hits in our local hit array
   // This will save CPU time for the isolation criterion 
   if (clean) hits.Add(sx);
  }
 }
 
 // Isolation cut
 // Only retain hits that have at least one hit of another OM within a certain
 // radius and within a certain time window
 Int_t nhits=hits.GetEntries();
 NcSignal* sx1=0;
 NcSignal* sx2=0;
 Float_t t1,t2;
 IceAOM* omx1=0;
 IceAOM* omx2=0;
 NcPosition r1;
 NcPosition r2;
 Float_t dt,dr;
 Int_t iso;
 for (Int_t jh1=0; jh1<nhits; jh1++)
 {
  sx1=(NcSignal*)hits.At(jh1);
  if (!sx1) continue;
 
  iso=1;
  for (Int_t jh2=0; jh2<nhits; jh2++)
  {
   sx2=(NcSignal*)hits.At(jh2);
   if (!sx2) continue;
 
   omx1=(IceAOM*)sx1->GetDevice();
   omx2=(IceAOM*)sx2->GetDevice();
   if (omx1==omx2) continue;
 
   t1=sx1->GetSignal("LE",7);
   t2=sx2->GetSignal("LE",7);
   dt=fabs(t2-t1);
   if (dt>fDtmaxA) continue;
 
   if (omx1 && omx2)
   {
    r1=omx1->GetPosition();
    r2=omx2->GetPosition();
    dr=r1.GetDistance(r2);
    if (dr>fRmaxA) continue;
    iso=0;
   }
  }
  if (iso) sx1->SetDead("LE");
 }   
}

void IceCleanHits::TWRDaq()
{
 
 // All Amanda OMs with a signal
 TObjArray* aoms=fEvt->GetDevices("IceAOM");
 if (!aoms) return;
 
 // Local OM array with bad/dead OMs (as indicated via IceCalibrate) discarded
 TObjArray oms;
 IceAOM* omx=0;
 for (Int_t i=0; i<aoms->GetEntries(); i++)
 {
  omx=(IceAOM*)aoms->At(i);
  if (!omx) continue;
  if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
  oms.Add(omx);
 }
 
 // Determination of the selected trigger time
 if (fTnamAT != "user") // Trigger time derived from trigger data
 {
  fTtimAT=fEvt->GetTriggerTime(fTnamAT);
 }
 else // Trigger time from median of hit times
 {
  if (fTtimAT<0)
  {
   fTtimAT=fEvt->GetTriggerTime(fTnamAT,&oms);
 
   // Make the value negative to stay with median of hit times selection
   // for all following invokations.
   if (fTtimAT>0) fTtimAT*=-1.;
   if (!fTtimAT) fTtimAT=-0.001;
  }
 }
 
 // Store the trigger time in the parameters device
 NcDevice* params=(NcDevice*)fEvt->GetDevice("IceCleanHits");
 if (params) params->SetSignal(fTtimAT,"TtimAT");
 
 // Local array with clean hits
 TObjArray hits;
 Int_t omid=0;
 Int_t clean=1;
 NcSignal* sx=0;
 Float_t adc,le,tot;
 Int_t readout;
 for (Int_t iom=0; iom<oms.GetEntries(); iom++)
 {
  omx=(IceAOM*)oms.At(iom);
  if (!omx) continue;
  omid=omx->GetUniqueID();
  readout=int(omx->GetSignal("READOUT"));
  // General readout setting in case info was missing 
  if (!readout)
  {
   readout=1;
   if (omid>=303) readout=2; // Optical OMs
  }
  for (Int_t ih=1; ih<=omx->GetNhits(); ih++)
  {
   sx=omx->GetHit(ih);
   if (!sx) continue;
   adc=sx->GetSignal("ADC",7);
   le=sx->GetSignal("LE",7);
   tot=sx->GetSignal("TOT",7);
 
   // Remove hits with an ADC value outside the range
   if (adc<fAdcminAT || adc>fAdcmaxAT)
   {
    sx->SetDead("ADC");
    clean=0;
   }
   // Remove hits with a TOT value outside the range
   // Note : Different ranges for electrical and optical modules
   if (readout==1) // Electrical OMs
   {
    if (tot<fTotminAET || tot>fTotmaxAET)
    {
     sx->SetDead("TOT");
     clean=0;
    }
   }
   else // Optical OMs
   {
    if (tot<fTotminAOT || tot>fTotmaxAOT)
    {
     sx->SetDead("TOT");
     clean=0;
    }
   }
   // Remove hits that are outside the trigger time window.
   // Since the trigger time was determined from uncalibrated LE's
   // (to include cable length effects) the uncalibrated LE of each
   // hit should be used here as well. 
   le=sx->GetSignal("LE",-7);
   if (fabs(le-fTtimAT)>fTwinAT)
   {
     sx->SetDead("LE");
     clean=0;
   }
   // Store only the current clean hits in our local hit array
   // This will save CPU time for the isolation criterion 
   if (clean) hits.Add(sx);
  }
 }
 
 // Isolation cut
 // Only retain hits that have at least one hit of another OM within a certain
 // radius and within a certain time window
 Int_t nhits=hits.GetEntries();
 NcSignal* sx1=0;
 NcSignal* sx2=0;
 Float_t t1,t2;
 IceAOM* omx1=0;
 IceAOM* omx2=0;
 NcPosition r1;
 NcPosition r2;
 Float_t dt,dr;
 Int_t iso;
 for (Int_t jh1=0; jh1<nhits; jh1++)
 {
  sx1=(NcSignal*)hits.At(jh1);
  if (!sx1) continue;
 
  iso=1;
  for (Int_t jh2=0; jh2<nhits; jh2++)
  {
   sx2=(NcSignal*)hits.At(jh2);
   if (!sx2) continue;
 
   omx1=(IceAOM*)sx1->GetDevice();
   omx2=(IceAOM*)sx2->GetDevice();
   if (omx1==omx2) continue;
 
   t1=sx1->GetSignal("LE",7);
   t2=sx2->GetSignal("LE",7);
   dt=fabs(t2-t1);
   if (dt>fDtmaxA) continue;
 
   if (omx1 && omx2)
   {
    r1=omx1->GetPosition();
    r2=omx2->GetPosition();
    dr=r1.GetDistance(r2);
    if (dr>fRmaxA) continue;
    iso=0;
   }
  }
  if (iso) sx1->SetDead("LE");
 }   
}

void IceCleanHits::IceTop()
{
 
 // All IceTop DOMs with a signal
 TObjArray* doms=fEvt->GetDevices("IceTDOM");
 if (!doms) return;
 
 TObjArray hits; // Local hit array
 IceGOM* omx=0;
 NcSignal* sx=0;
 for (Int_t i=0; i<doms->GetEntries(); i++)
 {
  omx=(IceGOM*)doms->At(i);
  if (!omx) continue;
  for (Int_t ih=1; ih<=omx->GetNhits(); ih++)
  {
   sx=omx->GetHit(ih);
   if (!sx) continue;
 
   if (omx->GetDeadValue("ADC")) sx->SetDead("ADC");
   if (omx->GetDeadValue("LE")) sx->SetDead("LE");
   if (omx->GetDeadValue("TOT")) sx->SetDead("TOT");
 
   hits.Add(sx);
  }
 }
 
 // ADC and TOT cleaning
 Amplitude(hits,fAdcminT,fAdcmaxT,fTotminT,fTotmaxT);
 
 // Isolated hit cleaning
 Isolation(hits,fDtmaxT,fRmaxT);
 
 // Determination of the selected trigger time
 TArrayF peaks;
 TArrayF* refp=0;
 if (fMultT) refp=&peaks;
 if (fTnamT != "user") // Trigger time derived from trigger data
 {
  fTtimT=fEvt->GetTriggerTime(fTnamT);
 }
 else // Trigger time from median of hit times
 {
  if (fTtimT<0)
  {
   fTtimT=fEvt->GetTriggerTime(fTnamT,&hits,fSlcT,refp);
 
   // Make the value negative to stay with median of hit times selection
   // for all following invokations.
   if (fTtimT>0) fTtimT*=-1.;
   if (!fTtimT) fTtimT=-0.001;
  }
 }
 
 // Store the trigger time in the parameters device
 NcDevice* params=(NcDevice*)fEvt->GetDevice("IceCleanHits");
 if (params) params->SetSignal(fTtimT,"TtimT");
 
 // Time window cleaning 
 TimeWindow(hits,fTtimT,fTwinT,refp);
}

void IceCleanHits::ICdoms()
{
 
 // All standard IceCube DOMs with a signal
 TObjArray* doms=fEvt->GetDevices("IceICDOM");
 if (!doms) return;
 
 TObjArray hits; // Local hit array
 IceGOM* omx=0;
 NcSignal* sx=0;
 for (Int_t i=0; i<doms->GetEntries(); i++)
 {
  omx=(IceGOM*)doms->At(i);
  if (!omx) continue;
  for (Int_t ih=1; ih<=omx->GetNhits(); ih++)
  {
   sx=omx->GetHit(ih);
   if (!sx) continue;
 
   if (omx->GetDeadValue("ADC")) sx->SetDead("ADC");
   if (omx->GetDeadValue("LE")) sx->SetDead("LE");
   if (omx->GetDeadValue("TOT")) sx->SetDead("TOT");
 
   hits.Add(sx);
  }
 }
 
 // ADC and TOT cleaning
 Amplitude(hits,fAdcminIC,fAdcmaxIC,fTotminIC,fTotmaxIC);
 
 // Isolated hit cleaning
 Isolation(hits,fDtmaxIC,fRmaxIC);
 
 // Determination of the selected trigger time
 TArrayF peaks;
 TArrayF* refp=0;
 if (fMultIC) refp=&peaks;
 if (fTnamIC != "user") // Trigger time derived from trigger data
 {
  fTtimIC=fEvt->GetTriggerTime(fTnamIC);
 }
 else // Trigger time from median of hit times
 {
  if (fTtimIC<0)
  {
   fTtimIC=fEvt->GetTriggerTime(fTnamIC,&hits,fSlcIC,refp);
 
   // Make the value negative to stay with median of hit times selection
   // for all following invokations.
   if (fTtimIC>0) fTtimIC*=-1.;
   if (!fTtimIC) fTtimIC=-0.001;
  }
 }
 
 // Store the trigger time in the parameters device
 NcDevice* params=(NcDevice*)fEvt->GetDevice("IceCleanHits");
 if (params) params->SetSignal(fTtimIC,"TtimIC");
 
 // Time window cleaning 
 TimeWindow(hits,fTtimIC,fTwinIC,refp);
}

void IceCleanHits::DCdoms()
{
 
 // All DeepCore DOMs with a signal
 TObjArray* doms=fEvt->GetDevices("IceDCDOM");
 if (!doms) return;
 
 TObjArray hits; // Local hit array
 IceGOM* omx=0;
 NcSignal* sx=0;
 for (Int_t i=0; i<doms->GetEntries(); i++)
 {
  omx=(IceGOM*)doms->At(i);
  if (!omx) continue;
  for (Int_t ih=1; ih<=omx->GetNhits(); ih++)
  {
   sx=omx->GetHit(ih);
   if (!sx) continue;
 
   if (omx->GetDeadValue("ADC")) sx->SetDead("ADC");
   if (omx->GetDeadValue("LE")) sx->SetDead("LE");
   if (omx->GetDeadValue("TOT")) sx->SetDead("TOT");
 
   hits.Add(sx);
  }
 }
 
 // ADC and TOT cleaning
 Amplitude(hits,fAdcminDC,fAdcmaxDC,fTotminDC,fTotmaxDC);
 
 // Isolated hit cleaning
 Isolation(hits,fDtmaxDC,fRmaxDC);
 
 // Determination of the selected trigger time
 TArrayF peaks;
 TArrayF* refp=0;
 if (fMultDC) refp=&peaks;
 if (fTnamDC != "user") // Trigger time derived from trigger data
 {
  fTtimDC=fEvt->GetTriggerTime(fTnamDC);
 }
 else // Trigger time from median of hit times
 {
  if (fTtimDC<0)
  {
   fTtimDC=fEvt->GetTriggerTime(fTnamDC,&hits,fSlcDC,refp);
 
   // Make the value negative to stay with median of hit times selection
   // for all following invokations.
   if (fTtimDC>0) fTtimDC*=-1.;
   if (!fTtimDC) fTtimDC=-0.001;
  }
 }
 
 // Store the trigger time in the parameters device
 NcDevice* params=(NcDevice*)fEvt->GetDevice("IceCleanHits");
 if (params) params->SetSignal(fTtimDC,"TtimDC");
 
 // Time window cleaning 
 TimeWindow(hits,fTtimDC,fTwinDC,refp);
}

void IceCleanHits::Amplitude(TObjArray& hits,Float_t adcmin,Float_t adcmax,Float_t totmin,Float_t totmax) const
{
 
 NcSignal* sx=0;
 Float_t adc,tot;
 for (Int_t ih=0; ih<hits.GetEntries(); ih++)
 {
  sx=(NcSignal*)hits.At(ih);
  if (!sx) continue;
  adc=sx->GetSignal("ADC",7);
  tot=sx->GetSignal("TOT",7);
 
  // Reject hits with an ADC value outside the range
  if (adc<adcmin || adc>adcmax) sx->SetDead("ADC");
 
  // Reject hits with a TOT value outside the range
  if (tot<totmin || tot>totmax) sx->SetDead("TOT");
 }
}

void IceCleanHits::Isolation(TObjArray& hits,Float_t dtmax,Float_t rmax) const
{
 
 NcSignal* sx1=0;
 NcSignal* sx2=0;
 Float_t t1,t2;
 IceGOM* omx1=0;
 IceGOM* omx2=0;
 NcPosition r1;
 NcPosition r2;
 Float_t dt,dr;
 Int_t iso;
 Int_t nhits=hits.GetEntries();
 for (Int_t jh1=0; jh1<nhits; jh1++)
 {
  sx1=(NcSignal*)hits.At(jh1);
  if (!sx1) continue;
 
  // No need to investigate further if LE was already rejected
  if (sx1->GetDeadValue("LE")) continue;
 
  iso=1;
 
  for (Int_t jh2=0; jh2<nhits; jh2++)
  {
   sx2=(NcSignal*)hits.At(jh2);
   if (!sx2) continue;
 
   if (sx2->GetDeadValue("ADC") || sx2->GetDeadValue("LE") || sx2->GetDeadValue("TOT")) continue;
 
   omx1=(IceGOM*)sx1->GetDevice();
   omx2=(IceGOM*)sx2->GetDevice();
   if (omx1==omx2) continue;
 
   t1=sx1->GetSignal("LE",7);
   t2=sx2->GetSignal("LE",7);
   dt=fabs(t2-t1);
   if (dt>dtmax) continue;
 
   if (omx1 && omx2)
   {
    r1=omx1->GetPosition();
    r2=omx2->GetPosition();
    dr=r1.GetDistance(r2);
    if (dr>rmax) continue;
    iso=0;
    break;
   }
  } // End of loop over hit2
 
  if (iso) sx1->SetDead("LE");
 } // End of loop over hit1   
}

void IceCleanHits::TimeWindow(TObjArray& hits,Float_t ttrig,Float_t twin,TArrayF* peaks) const
{
 
 NcSignal* sx=0;
 Float_t le;
 Int_t npeaks=0;
 if (peaks) npeaks=peaks->GetSize();
 if (npeaks<1) npeaks=1;
 for (Int_t ih=0; ih<hits.GetEntries(); ih++)
 {
  sx=(NcSignal*)hits.At(ih);
  if (!sx) continue;
 
  // No need to investigate further if LE was already rejected
  if (sx->GetDeadValue("LE")) continue;
 
  le=sx->GetSignal("LE",7);
 
  // Retain only hits that are inside a trigger time window.
  sx->SetDead("LE");
  for (Int_t ipeak=0; ipeak<npeaks; ipeak++)
  {
   if (peaks) ttrig=peaks->At(ipeak);
   if (fabs(le-fabs(ttrig))<=twin)
   {
    sx->SetAlive("LE");
    break;
   }
  }
 }
}