RALICE/icepack/IceXtalk.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 IceXtalk
  20 // TTask derived class to perform cross talk hit correction.
  21 //
  22 // This task takes the current event in memory and uses the attached
  23 // OM database to identify pairs of OMs which might induce cross talk.
  24 // For each particular transmitter and receiver pair within the event
  25 // a probability for cross talk induction is determined on basis of
  26 // the hit data and the cross talk probability functions from the OM database.
  27 // If this probability is above a certain threshold "pmin" and the difference
  28 // in LE between transmitter and receiver signal is within the bounds as
  29 // given by the Xtalk calibration data, the occurrence of a cross talk signal
  30 // of "pe" photo-electrons in the receiver is assumed and this signal is
  31 // lateron subtracted from the receiver ADC signal.
  32 // The subtraction of the Xtalk signal from the various hits is delayed
  33 // until the very end, since receiver OMs may also act as transmitter OMs
  34 // for other pair combinations.
  35 // In case a transmitter OM was flagged as a bad module, the check on the
  36 // LE time difference will not be performed and a cross talk induction
  37 // decision will be solely based on the corresponding probability, since
  38 // the latter only depends on the un-calibrated signal amplitude of the
  39 // transmitter hit.
  40 // So, hits which have induced cross talk are not completely removed
  41 // but their signal is corrected for the cross talk.
  42 // This will prevent severe losses when studying UHE events, where cross
  43 // talk effects are expected to have negligible effects on the observed
  44 // large module signals.
  45 //
  46 // Note : The amount with which the ADC value was corrected is stored
  47 //        in the signal as an ADC offset. This will allow easy investigation
  48 //        of Xtalk corrected signals and also enables successive
  49 //        applications of this Xtalk processor to investigate (systematic)
  50 //        effects of various criteria.
  51 //        Example : In case an amount of 1 pe was subtracted from the ADC
  52 //                  value, the ADC offset will be set to 1.
  53 //
  54 // The values of "pmin" and "pe" can be set by the user via the
  55 // SetMinProb() and SetXtalkPE() memberfunctions.
  56 // The default values are pmin=0.5 and pe=1.
  57 //
  58 // In case this processor is followed by an ADC threshold hit cleaning
  59 // procedure, hits which only contained cross talk can be efficiently
  60 // skipped or removed from the event.
  61 // In case one would like to always remove a hit which containes induced
  62 // cross talk, one could set the "pe" parameter to some very large value.
  63 // This will result in cross talk induced hits to get large negative ADC
  64 // signals and can as such easily be skipped/removed afterwards.
  65 //
  66 // Note : This processor only works properly on Time and ADC calibrated data.
  67 //        In case no OM database has been specified for this processor,
  68 //        no cross talk hit correction will be performed.
  69 //
  70 //--- Author: Nick van Eijndhoven 11-aug-2005 Utrecht University
  71 //- Modified: NvE $Date$ Utrecht University
  72 ///////////////////////////////////////////////////////////////////////////
  73
  74 #include "IceXtalk.h"
  75 #include "Riostream.h"
  76
  77 ClassImp(IceXtalk) // Class implementation to enable ROOT I/O
  78
  79 IceXtalk::IceXtalk(const char* name,const char* title) : TTask(name,title)
  80 {
  81 // Default constructor.
  82  fCalfile=0;
  83  fOmdb=0;
  84  fPmin=0.5;
  85  fPe=1;
  86 }
  87 ///////////////////////////////////////////////////////////////////////////
  88 IceXtalk::~IceXtalk()
  89 {
  90 // Default destructor.
  91  if (fCalfile)
  92  {
  93   delete fCalfile;
  94   fCalfile=0;
  95  }
  96 }
  97 ///////////////////////////////////////////////////////////////////////////
  98 void IceXtalk::SetOMdbase(AliObjMatrix* omdb)
  99 {
 100 // Set the pointer to the OM database.
 101 // Note : this will overrule a previously attached database.
 102  fOmdb=omdb;
 103 }
 104 ///////////////////////////////////////////////////////////////////////////
 105 void IceXtalk::SetCalibFile(TString name)
 106 {
 107 // Set the calibration ROOT file as created with IceCal2Root.
 108 // Note : this will overrule a previously attached database.
 109  fCalfile=new TFile(name.Data());
 110  if (fCalfile)
 111  {
 112   fOmdb=(AliObjMatrix*)fCalfile->Get("Cal-OMDBASE");
 113  }
 114 }
 115 ///////////////////////////////////////////////////////////////////////////
 116 void IceXtalk::SetMinProb(Float_t pmin)
 117 {
 118 // Set the minimal probability for cross talk induction.
 119  fPmin=pmin;
 120 }
 121 ///////////////////////////////////////////////////////////////////////////
 122 void IceXtalk::SetXtalkPE(Float_t pe)
 123 {
 124 // Set the nominal Xtalk signal in photo-electron equivalent.
 125  fPe=pe;
 126 }
 127 ///////////////////////////////////////////////////////////////////////////
 128 void IceXtalk::Exec(Option_t* opt)
 129 {
 130 // Implementation of cross talk hit correction.
 131
 132  if (!fOmdb) return;
 133
 134  TString name=opt;
 135  AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
 136
 137  if (!parent) return;
 138
 139  IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent");
 140  if (!evt) return;
 141
 142  // All Amanda OMs with a signal
 143  TObjArray* mods=evt->GetDevices("IceAOM");
 144
 145  Int_t nmods=mods->GetEntries();
 146  if (!nmods) return;
 147
 148  TObjArray xhits; // Array with pointers to Xtalk hits to be corrected
 149
 150  IceAOM* omt=0; // Transmitter OM
 151  IceAOM* omr=0; // Receiver OM
 152  TF1* fxtalk=0; // The Xtalk probability function
 153  Int_t idtrans=0;
 154  Int_t idrec=0;
 155  Float_t adct=0,let=0;
 156  Float_t p=0;
 157  Float_t adcr=0,ler=0;
 158  Int_t ibad=0;
 159  Float_t cpar=0,bpar=0,dlemin=0,dlemax=0,test=0;
 160  Float_t dle=0;
 161  Float_t sigcor=0;
 162  for (Int_t imod=0; imod<nmods; imod++)
 163  {
 164   omt=(IceAOM*)mods->At(imod);
 165   if (!omt) continue;
 166   idtrans=omt->GetUniqueID();
 167
 168   // Also take bad transmitter modules into account
 169   ibad=omt->GetDeadValue("ADC");
 170   if (ibad) omt->SetAlive("ADC");
 171
 172   // Check for corresponding receiver modules
 173   for (Int_t jmod=0; jmod<nmods; jmod++)
 174   {
 175    // No Xtalk from a module to itself
 176    if (jmod==imod) continue;
 177
 178    omr=(IceAOM*)mods->At(jmod);
 179    if (!omr) continue;
 180    idrec=omr->GetUniqueID();
 181
 182    fxtalk=(TF1*)fOmdb->GetObject(idtrans,idrec+1);
 183    if (!fxtalk) continue;
 184
 185    // Determine Xtalk probability for each transmitter hit
 186    for (Int_t ithit=1; ithit<=omt->GetNhits(); ithit++)
 187    {
 188     AliSignal* st=omt->GetHit(ithit);
 189     if (!st) continue;
 190
 191     // Eliminate possible previous Xtalk correction
 192     sigcor=st->GetOffset("ADC");
 193     st->AddSignal(sigcor,"ADC");
 194     st->ResetOffset("ADC");
 195     adct=st->GetSignal("ADC",-4); // Get uncalibrated amplitude
 196     dlemin=fxtalk->GetParameter("dLE-min");
 197     dlemax=fxtalk->GetParameter("dLE-max");
 198     // Protection against overflow in Xtalk prob. function
 199     cpar=fxtalk->GetParameter("C");
 200     bpar=fxtalk->GetParameter("B");
 201     test=(adct-cpar)/bpar;
 202     if (test>100)
 203     {
 204      p=1;
 205     }
 206     else if (test<-100)
 207     {
 208      p=0;
 209     }
 210     else
 211     {
 212      p=fxtalk->Eval(adct);
 213     }
 214     if (p<fPmin) continue;
 215
 216     // Xtalk flagging for each receiver hit
 217     for (Int_t irhit=1; irhit<=omr->GetNhits(); irhit++)
 218     {
 219      AliSignal* sr=omr->GetHit(irhit);
 220      if (!sr) continue;
 221
 222      // Check calibrated LE time differences
 223      if (!ibad)
 224      {
 225       let=st->GetSignal("LE",4);
 226       ler=sr->GetSignal("LE",4);
 227       dle=ler-let;
 228       if (dle<dlemin || dle>dlemax) continue;
 229      }
 230      // Register this receiver hit to be corrected
 231      xhits.Add(sr);
 232     }
 233    } // end of transmitter hit loop
 234   } // end of receiver OM loop
 235
 236   // Restore the original transmitter dead flag value
 237   if (ibad) omt->SetDead("ADC");
 238  } // end of transmitter OM loop
 239
 240  // Perform the Xtalk signal correction to the registered receiver hits
 241  // The correction value is stored in the signal data as an offset value
 242  for (Int_t ix=0; ix<xhits.GetEntries(); ix++)
 243  {
 244   AliSignal* sx=(AliSignal*)xhits.At(ix);
 245   if (!sx) continue;
 246
 247   // Eliminate possible previous Xtalk correction
 248   sigcor=sx->GetOffset("ADC");
 249   adcr=sx->GetSignal("ADC")+sigcor;
 250   sigcor=fPe;
 251   // If stored ADC data is un-calibrated, convert fPe to raw ADC value
 252   AliSignal* parent=(AliSignal*)sx->GetDevice();
 253   if (parent)
 254   {
 255    if (parent->GetCalFunction("ADC"));
 256    {
 257     idrec=parent->GetUniqueID();
 258     omr=(IceAOM*)fOmdb->GetObject(idrec,1);
 259     TF1* fdecal=0;
 260     if (omr) fdecal=omr->GetDecalFunction("ADC");
 261     if (fdecal) sigcor=fdecal->Eval(fPe);
 262    }
 263   }
 264   adcr=adcr-sigcor;
 265   sx->SetSignal(adcr,"ADC");
 266   sx->SetOffset(sigcor,"ADC");
 267  }
 268 }
 269 ///////////////////////////////////////////////////////////////////////////