RALICE/icepack/IceMakeHits.cxx

   1 /*******************************************************************************
   2  * Copyright(c) 2003, IceCube Experiment at the South Pole. All rights reserved.
   3  *
   4  * Author: The IceCube RALICE-based Offline Project.
   5  * Contributors are mentioned in the code where appropriate.
   6  *
   7  * Permission to use, copy, modify and distribute this software and its
   8  * documentation strictly for non-commercial purposes is hereby granted
   9  * without fee, provided that the above copyright notice appears in all
  10  * copies and that both the copyright notice and this permission notice
  11  * appear in the supporting documentation.
  12  * The authors make no claims about the suitability of this software for
  13  * any purpose. It is provided "as is" without express or implied warranty.
  14  *******************************************************************************/
  15
  16 // $Id$
  17
  18 ///////////////////////////////////////////////////////////////////////////
  19 // Class IceMakeHits
  20 // TTask derived class to perform hit extraction from waveforms.
  21 //
  22 // The code in this processor is based on the algorithms as developed by
  23 // Nick van Eijndhoven and Garmt de Vries-Uiterweerd (Utrecht University, The Netherlands).
  24 //
  25 // Procedure applied for Amanda TWR data :
  26 // ---------------------------------------
  27 //
  28 // 1) The waveform is fed to a TSpectrum object, and the peak locations
  29 //    are determined with the TSpectrum::Search() function.
  30 //
  31 // 2) The waveform is divided into regions corresponding to the peaks found by
  32 //    TSpectrum. The region boundary between two peaks is at the location of
  33 //    the minimum between the two peaks.
  34 //
  35 // 3) For each region the "effective baseline" (used in the
  36 //    evaluation of the leading edge value) is determined as :
  37 //    effective baseline = fBasefracXXX * value at lower region boundary.
  38 //    This takes into account the effect from the previous pulse.
  39 //    For the first pulse, the effective baseline is equal to the overall
  40 //    baseline.
  41 //
  42 // 4) For each region, the point of steepest rise between the lower region
  43 //    boundary and the peak location is determined. The tangent at this point
  44 //    is extrapolated to the effective baseline. The point of intersection yields the
  45 //    leading edge.
  46 //
  47 // 5) For each region the range of charge integration is determined as :
  48 //    - Start of integration at the lower region boundary or at the leading edge,
  49 //      whichever comes last;
  50 //    - End of integration at the upper region boundary or at the point where the
  51 //      signal drops below the overall baseline, whichever comes first.
  52 //
  53 // 6) For each region the integrated charge is determined as :
  54 //    Sum over bins in integration range of (value in bin - overall baseline).
  55 //
  56 // 7) For each pulse the quality is evaluated by requiring that :
  57 //    peak location - lower region boundary > lower region boundary - leading edge.
  58 //    For a too shallow steepest rise, the leading edge value is unreliable, in
  59 //    which case the pulse is merged with the previous pulse.
  60 //
  61 // 8) Each pulse is checked for saturation and discarded if necessary.
  62 //
  63 // 9) TSpectrum needs a minimum number of bins for its Search function, otherwise
  64 //    the clipping window is too large, which causes an error. If a waveform does
  65 //    not contain enough bins, the following alternative approach is used :
  66 //    - A loop over all bins is performed.
  67 //    - As soon as the signal exceeds a given threshold, a pulse is started.
  68 //    - The pulse ends when the signal drops below the threshold again.
  69 //    - While looping, the charge is integrated for each pulse.
  70 //
  71 // The defaults of the various parameters can be changed by the corresponding
  72 // Set memberfunctions.
  73 //
  74 // Information about the actual parameter settings can be found in the event
  75 // structure itself via the device named "IceMakeHits".
  76 //
  77 //--- Author: Nick van Eijndhoven and Garmt de Vries-Uiterweerd 15-jan-2007 Utrecht University
  78 //- Modified: NvE $Date$ Utrecht University
  79 ///////////////////////////////////////////////////////////////////////////
  80
  81 #include "IceMakeHits.h"
  82 #include "Riostream.h"
  83
  84 ClassImp(IceMakeHits) // Class implementation to enable ROOT I/O
  85
  86 IceMakeHits::IceMakeHits(const char* name,const char* title) : TTask(name,title)
  87 {
  88 // Default constructor.
  89  fEvt=0;
  90  fDaq=0;
  91  fBasefracA=0.5;
  92  fSigmaA=1.5;
  93  fMaxPeaksA=10;
  94  fMinPulseHeightA=50;
  95  fThresholdA=0.2;
  96 }
  97 ///////////////////////////////////////////////////////////////////////////
  98 IceMakeHits::~IceMakeHits()
  99 {
 100 // Default destructor.
 101 }
 102 ///////////////////////////////////////////////////////////////////////////
 103 void IceMakeHits::SetBasefracA(Float_t val)
 104 {
 105 // Set baseline fractional update for Amanda TWR extraction.
 106 // The default as set in the constructor of this class is 0.5.
 107  fBasefracA=val;
 108 }
 109 ///////////////////////////////////////////////////////////////////////////
 110 void IceMakeHits::SetSigmaA(Float_t val)
 111 {
 112 // Set clipping window width for Amanda TWR extraction.
 113 // The default as set in the constructor of this class is 1.5.
 114  fSigmaA=val;
 115 }
 116 ///////////////////////////////////////////////////////////////////////////
 117 void IceMakeHits::SetMaxPeaksA(Int_t val)
 118 {
 119 // Set maximum number of peaks in a waveform for Amanda TWR extraction.
 120 // The default as set in the constructor of this class is 10.
 121  fMaxPeaksA=val;
 122 }
 123 ///////////////////////////////////////////////////////////////////////////
 124 void IceMakeHits::SetMinPulseHeightA(Float_t val)
 125 {
 126 // Set minimum required pulse height for Amanda TWR extraction.
 127 // This is used only for narrow pulses that cannot be handled with TSpectrum.
 128 // The default as set in the constructor of this class is 50.
 129  fMinPulseHeightA=val;
 130 }
 131 ///////////////////////////////////////////////////////////////////////////
 132 void IceMakeHits::SetThresholdA(Float_t val)
 133 {
 134 // Set threshold for use in analysis of narrow pulses for Amanda TWR extraction.
 135 // A peak is assumed to start when the signal rises above threshold*maxval,
 136 // where maxval is the maximum value found in the waveform.
 137 // The default as set in the constructor of this class is 0.2.
 138  fThresholdA=val;
 139 }
 140 ///////////////////////////////////////////////////////////////////////////
 141 void IceMakeHits::Exec(Option_t* opt)
 142 {
 143 // Implementation of the hit cleaning procedures.
 144
 145  TString name=opt;
 146  AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
 147
 148  if (!parent) return;
 149
 150  fEvt=(IceEvent*)parent->GetObject("IceEvent");
 151  if (!fEvt) return;
 152
 153  fDaq=(AliDevice*)fEvt->GetDevice("Daq");
 154  if (!fDaq) return;
 155
 156  // Storage of the used parameters in the IceMakeHits device
 157  AliSignal params;
 158  params.SetNameTitle("IceMakeHits","IceMakeHits processor parameters");
 159  params.SetSlotName("BasefracA",1);
 160  params.SetSlotName("SigmaA",2);
 161  params.SetSlotName("MaxPeaksA",3);
 162  params.SetSignal(fBasefracA,1);
 163  params.SetSignal(fSigmaA,2);
 164  params.SetSignal(fMaxPeaksA,3);
 165
 166  fEvt->AddDevice(params);
 167
 168  Amanda();
 169  InIce();
 170  IceTop();
 171 }
 172 ///////////////////////////////////////////////////////////////////////////
 173 void IceMakeHits::Amanda()
 174 {
 175 // Hit extraction from the Amanda TWR data.
 176
 177  // Arrays for storing info
 178  Float_t* baseline=new Float_t[fMaxPeaksA];
 179  Int_t* lowend=new Int_t[fMaxPeaksA];
 180  Int_t* upend=new Int_t[fMaxPeaksA];
 181  Int_t* startcharge=new Int_t[fMaxPeaksA];
 182  Int_t* stopcharge=new Int_t[fMaxPeaksA];
 183  Int_t* status=new Int_t[fMaxPeaksA]; // 0=OK, 1=rejected, 2=saturation
 184  Float_t* leadingedge=new Float_t[fMaxPeaksA];
 185  Float_t* charge=new Float_t[fMaxPeaksA];
 186  Float_t* tot=new Float_t[fMaxPeaksA];
 187
 188  // Some objects and variables we will need
 189  TH1F foo, diff;
 190  TSpectrum* spec=new TSpectrum(fMaxPeaksA);
 191  Int_t nrIterations=(Int_t)(7*fSigmaA+0.5); // Number of iterations used in TSpectrum::SearchHighRes()
 192  Int_t npeaks=0, ibin=0, lookforsteepestuntilbin=0, steep=0;
 193  Float_t maxval=0, rise=0, rc=0, yyy=0;
 194  Bool_t pulsegoingon=false;
 195  Int_t* index=new Int_t[fMaxPeaksA];
 196
 197  // Update the DAQ device data for the nature of the hit info
 198  fDaq->AddNamedSlot("TWR");
 199  fDaq->SetSignal(1,"TWR");
 200  Int_t idx=fDaq->GetSlotIndex("Muon");
 201  if (idx) fDaq->SetSignal(0,idx);
 202
 203  // All Amanda OMs with a signal
 204  TObjArray* aoms=fEvt->GetDevices("IceAOM");
 205
 206  // OM, waveform and hit
 207  IceAOM* omx=0;
 208  TH1F* wf=0;
 209  AliSignal hit;
 210  hit.SetSlotName("ADC",1);
 211  hit.SetSlotName("LE",2);
 212  hit.SetSlotName("TOT",3);
 213
 214  // Loop over all fired OMs and extract the hit info
 215  for (Int_t iom=0; iom<aoms->GetEntries(); iom++)
 216  {
 217   omx=(IceAOM*)aoms->At(iom);
 218   if (!omx) continue;
 219   // Remove all existing hits of this OM
 220   omx->RemoveHits();
 221
 222   // Should we skip OMs that we know from the dbase to have problems ?
 223 ////  if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
 224
 225   // Investigate all waveforms for this OM
 226   for (Int_t iwf=1; iwf<=omx->GetNwaveforms(); iwf++)
 227   {
 228    wf=omx->GetWaveform(iwf);
 229    if (!wf) continue;
 230
 231    maxval=wf->GetMaximum();
 232
 233    // Check if clipping window is not too large
 234    if(wf->GetNbinsX() > 2*nrIterations+1)
 235    {
 236     // Find peaks with TSpectrum
 237     npeaks=spec->Search(wf,fSigmaA,"goff");
 238
 239     // Discard waveform if no or too many peaks found
 240     if(npeaks<1 || npeaks>fMaxPeaksA) continue;
 241
 242     // Get differential of WF
 243     diff=*wf;
 244     for(ibin=2;ibin<diff.GetNbinsX();ibin++)
 245     {
 246      diff.SetBinContent(ibin,wf->GetBinContent(ibin)-wf->GetBinContent(ibin-1));
 247     }
 248     diff.SetBinContent(1,0);
 249
 250     // Set baseline and lower end for first peak,
 251     baseline[0]=0;
 252     lowend[0]=1;
 253
 254     // Sort peaks in time
 255     TMath::Sort(npeaks,spec->GetPositionX(),index,false);
 256     // For each of the peaks,
 257     for(Int_t ipeak=0; ipeak<npeaks; ipeak++)
 258     {
 259      // Find baseline and region around peak
 260      foo=*wf;
 261      // (Second and later peaks: lower edge = upper edge previous peak,
 262      // baseline is average of previous baseline and minimum value between two
 263      // peaks)
 264      if(ipeak>0)
 265      {
 266       lowend[ipeak]=upend[ipeak-1]+1;
 267       baseline[ipeak]=fBasefracA*foo.GetBinContent(lowend[ipeak]);
 268      }
 269      // (Upper edge range is minimum between this and next peak)
 270      if(ipeak<npeaks-1)
 271      {
 272       foo.SetAxisRange(spec->GetPositionX()[index[ipeak]],spec->GetPositionX()[index[ipeak+1]]);
 273       upend[ipeak]=foo.GetMinimumBin();
 274      }
 275      // (Last peak: upper edge is end of histo)
 276      else
 277      {
 278       upend[ipeak]=wf->GetNbinsX();
 279      }
 280
 281      // Find steepest rise
 282      lookforsteepestuntilbin=wf->FindBin(spec->GetPositionX()[index[ipeak]]);
 283      foo=diff;
 284      foo.SetAxisRange(wf->GetBinLowEdge(lowend[ipeak]),wf->GetBinLowEdge(lookforsteepestuntilbin+1));
 285      do
 286      {
 287       steep=foo.GetMaximumBin();
 288       rise=foo.GetBinContent(steep);
 289       if(rise==1e-9) break;
 290       foo.SetBinContent(steep,-1e9);
 291      } while(wf->GetBinContent(steep)<baseline[ipeak]);
 292
 293      // Extrapolate tangent to find leading edge
 294      yyy=wf->GetBinContent(steep)-baseline[ipeak];
 295      rc=rise/foo.GetBinWidth(steep);
 296      if(rc>0) leadingedge[ipeak]=wf->GetBinCenter(steep)-yyy/rc; else leadingedge[ipeak]=0;
 297
 298      // Determine peak status
 299      status[ipeak]=0;
 300      // Check for saturation
 301      if(rc<0.1 && wf->GetBinContent(wf->FindBin(spec->GetPositionX()[index[ipeak]])) == maxval)
 302      {
 303       status[ipeak]=2;
 304      }
 305      // Check quality: LE should not be too far below lower edge
 306      // Otherwise, ignore this peak and set baseline back to what it was
 307      else if(wf->GetBinLowEdge(lowend[ipeak]) - leadingedge[ipeak] > spec->GetPositionX()[index[ipeak]] - wf->GetBinLowEdge(lowend[ipeak]))
 308      {
 309       status[ipeak]=1;
 310       if(ipeak>0) baseline[ipeak]=baseline[ipeak-1];
 311      }
 312
 313      // Start charge integration at LE, or at lower edge of range
 314      startcharge[ipeak]=wf->FindBin(leadingedge[ipeak]);
 315      if(lowend[ipeak]>startcharge[ipeak]) startcharge[ipeak]=lowend[ipeak];
 316
 317      // Integrate charge until pulse drop below baseline, or else until edge of range
 318      stopcharge[ipeak]=upend[ipeak];
 319      for(ibin=wf->FindBin(spec->GetPositionX()[index[ipeak]]); ibin<=upend[ipeak]; ibin++)
 320      {
 321       if(wf->GetBinContent(ibin)<0)
 322       {
 323        stopcharge[ipeak]=ibin-1;
 324        break;
 325       }
 326      }
 327
 328      // Determine time over threshold
 329      tot[ipeak]=wf->GetBinLowEdge(stopcharge[ipeak]+1)-wf->GetBinLowEdge(startcharge[ipeak]);
 330
 331      // Determine charge
 332      charge[ipeak]=0;
 333      for(ibin=startcharge[ipeak]; ibin<=stopcharge[ipeak]; ibin++)
 334      {
 335       charge[ipeak]+=wf->GetBinContent(ibin);
 336      }
 337
 338     } // end loop over peaks
 339
 340     // Check all peaks, from latest to earliest
 341     for(int ipeak=npeaks-1; ipeak>=0; ipeak--)
 342     {
 343
 344      // If this peak was rejected, add charge and TOT to previous peak (if there is one)
 345      if(status[ipeak]==1 && ipeak>0)
 346      {
 347       charge[ipeak-1]+=charge[ipeak];
 348       charge[ipeak]=0;
 349       tot[ipeak-1]+=tot[ipeak];
 350       tot[ipeak]=0;
 351      }
 352
 353      // If this peak is OK, add hit info
 354      if(status[ipeak]==0)
 355      {
 356       hit.Reset();
 357       hit.SetSignal(charge[ipeak],"ADC");
 358       hit.SetSignal(leadingedge[ipeak],"LE");
 359       hit.SetSignal(tot[ipeak],"TOT");
 360       omx->AddHit(hit);
 361      }
 362
 363     } // end loop over peaks
 364
 365    // If number of bins too small, use different method
 366    }
 367    else
 368    {
 369     // If maximum value high enough to suspect presence of peak,
 370     if(maxval>fMinPulseHeightA)
 371     {
 372      // Loop over bins
 373      pulsegoingon=false;
 374      npeaks=0;
 375      for(ibin=1; ibin<=wf->GetNbinsX(); ibin++)
 376      {
 377       // If bin content above threshold, start pulse
 378       if(wf->GetBinContent(ibin)>fThresholdA*maxval){
 379        if(!pulsegoingon)
 380        {
 381         // Pulse starts here
 382         pulsegoingon=true;
 383         leadingedge[npeaks]=wf->GetBinLowEdge(ibin);
 384         charge[npeaks]=wf->GetBinContent(ibin);
 385        }
 386        else
 387        {
 388         // Pulse continues
 389         charge[npeaks]+=wf->GetBinContent(ibin);
 390        }
 391       }
 392       else
 393       {
 394        if(pulsegoingon)
 395        {
 396         // Pulse ends here
 397         tot[npeaks]=wf->GetBinLowEdge(ibin)-leadingedge[npeaks];
 398
 399         // Store pulse information
 400         hit.Reset();
 401         hit.SetSignal(charge[npeaks],"ADC");
 402         hit.SetSignal(leadingedge[npeaks],"LE");
 403         hit.SetSignal(tot[npeaks],"TOT");
 404         omx->AddHit(hit);
 405
 406         // Get ready for next pulse
 407         pulsegoingon=false;
 408         npeaks++;
 409        }
 410       }
 411
 412      } // End of loop over bins
 413     }
 414
 415    } // End of alternative method for narrow pulses
 416
 417   } // End of WF loop
 418  } // End of OM loop
 419
 420  // Clean up
 421  delete[] baseline;
 422  delete[] lowend;
 423  delete[] upend;
 424  delete[] startcharge;
 425  delete[] stopcharge;
 426  delete[] status;
 427  delete[] leadingedge;
 428  delete[] charge;
 429  delete[] tot;
 430  delete[] index;
 431 }
 432 ///////////////////////////////////////////////////////////////////////////
 433 void IceMakeHits::InIce()
 434 {
 435 // Hit extraction from IceCube InIce ATWD data.
 436 }
 437 ///////////////////////////////////////////////////////////////////////////
 438 void IceMakeHits::IceTop()
 439 {
 440 // Hit extraction from IceTop ATWD data.
 441 }
 442 ///////////////////////////////////////////////////////////////////////////