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