12-aug-2005 NvE TTask derived class IceXtalk introduced for cross talk correction.
[u/mrichter/AliRoot.git] / 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 // The values of "pmin" and "pe" can be set by the user via the
47 // SetMinProb() and SetXtalkPE() memberfunctions.
48 // The default values are pmin=0.5 and pe=1.
49 //
50 // In case this processor is followed by an ADC threshold hit cleaning
51 // procedure, hits which only contained cross talk can be efficiently
52 // skipped or removed from the event.
53 // In case one would like to always remove a hit which containes induced
54 // cross talk, one could set the "pe" parameter to some very large value.
55 // This will result in cross talk induced hits to get large negative ADC
56 // signals and can as such easily be skipped/removed afterwards. 
57 //
58 // Note : This processor only works properly on Time and ADC calibrated data.
59 //
60 //--- Author: Nick van Eijndhoven 11-aug-2005 Utrecht University
61 //- Modified: NvE $Date$ Utrecht University
62 ///////////////////////////////////////////////////////////////////////////
63  
64 #include "IceXtalk.h"
65 #include "Riostream.h"
66
67 ClassImp(IceXtalk) // Class implementation to enable ROOT I/O
68
69 IceXtalk::IceXtalk(const char* name,const char* title) : TTask(name,title)
70 {
71 // Default constructor.
72  fCalfile=0;
73  fOmdb=0;
74  fPmin=0.5;
75  fPe=1;
76 }
77 ///////////////////////////////////////////////////////////////////////////
78 IceXtalk::~IceXtalk()
79 {
80 // Default destructor.
81  if (fCalfile)
82  {
83   delete fCalfile;
84   fCalfile=0;
85  }
86 }
87 ///////////////////////////////////////////////////////////////////////////
88 void IceXtalk::SetOMdbase(AliObjMatrix* omdb)
89 {
90 // Set the pointer to the OM database.
91 // Note : this will overrule a previously attached database. 
92  fOmdb=omdb;
93 }
94 ///////////////////////////////////////////////////////////////////////////
95 void IceXtalk::SetCalibFile(TString name)
96 {
97 // Set the calibration ROOT file as created with IceCal2Root.
98 // Note : this will overrule a previously attached database. 
99  fCalfile=new TFile(name.Data());
100  if (fCalfile)
101  {
102   AliObjMatrix* fOmdb=(AliObjMatrix*)fCalfile->Get("Cal-OMDBASE");
103  }
104 }
105 ///////////////////////////////////////////////////////////////////////////
106 void IceXtalk::SetMinProb(Float_t pmin)
107 {
108 // Set the minimal probability for cross talk induction.
109  fPmin=pmin;
110 }
111 ///////////////////////////////////////////////////////////////////////////
112 void IceXtalk::SetXtalkPE(Float_t pe)
113 {
114 // Set the nominal Xtalk signal in photo-electron equivalent.
115  fPe=pe;
116 }
117 ///////////////////////////////////////////////////////////////////////////
118 void IceXtalk::Exec(Option_t* opt)
119 {
120 // Implementation of cross talk hit correction.
121
122  TString name=opt;
123  AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
124
125  if (!parent) return;
126
127  IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent");
128  if (!evt) return;
129
130  TObjArray xhits; // Array with pointers to Xtalk hits to be corrected
131
132  IceAOM* omt=0; // Transmitter OM
133  IceAOM* omr=0; // Receiver OM
134  TF1* fxtalk=0; // The Xtalk probability function
135  Int_t idtrans=0;
136  Int_t idrec=0;
137  Float_t adct=0,let=0;
138  Float_t p=0;
139  Float_t adcr=0,ler=0;
140  Int_t ibad=0;
141  Float_t cpar=0,bpar=0,dlemin=0,dlemax=0,test=0;
142  Float_t dle=0;
143  Float_t sigcor=0;
144  Int_t ndev=evt->GetNdevices();
145  for (Int_t idev=1; idev<=ndev; idev++)
146  {
147   omt=(IceAOM*)evt->GetDevice(idev);
148   if (!omt) continue;
149   idtrans=omt->GetUniqueID();
150
151   // Also take bad transmitter modules into account
152   ibad=omt->GetDeadValue("ADC");
153   if (ibad) omt->SetAlive("ADC");
154
155   // Check for corresponding receiver modules  
156   for (Int_t jdev=1; jdev<=ndev; jdev++)
157   {
158    // No Xtalk from a module to itself 
159    if (jdev==idev) continue;
160
161    omr=(IceAOM*)evt->GetDevice(jdev);
162    if (!omr) continue;
163    idrec=omr->GetUniqueID();
164
165    fxtalk=(TF1*)fOmdb->GetObject(idtrans,idrec+1);
166    if (!fxtalk) continue;
167
168    // Determine Xtalk probability for each transmitter hit
169    for (Int_t ithit=1; ithit<=omt->GetNhits(); ithit++)
170    {
171     AliSignal* st=omt->GetHit(ithit);
172     if (!st) continue;
173
174     adct=st->GetSignal("ADC",-4); // Get uncalibrated amplitude
175     dlemin=fxtalk->GetParameter("dLE-min");
176     dlemax=fxtalk->GetParameter("dLE-max");
177     // Protection against overflow in Xtalk prob. function
178     cpar=fxtalk->GetParameter("C");
179     bpar=fxtalk->GetParameter("B");
180     test=(adct-cpar)/bpar;
181     if (test>100)
182     {
183      p=1;
184     }
185     else if (test<-100)
186     {
187      p=0;
188     }
189     else
190     {
191      p=fxtalk->Eval(adct);
192     }
193     if (p<fPmin) continue;
194
195     // Xtalk flagging for each receiver hit
196     for (Int_t irhit=1; irhit<=omr->GetNhits(); irhit++)
197     {
198      AliSignal* sr=omr->GetHit(irhit);
199      if (!sr) continue;
200
201      // Check calibrated LE time differences
202      if (!ibad)
203      {
204       let=st->GetSignal("LE",4);
205       ler=sr->GetSignal("LE",4);
206       dle=ler-let;
207       if (dle<dlemin || dle>dlemax) continue;
208      }
209      // Register this receiver hit to be corrected
210      xhits.Add(sr); 
211     }
212    } // end of transmitter hit loop
213   } // end of receiver OM loop
214
215   // Restore the original transmitter dead flag value
216   if (ibad) omt->SetDead("ADC");
217  } // end of transmitter OM loop
218  
219  // Perform the Xtalk signal correction to the registered receiver hits
220  for (Int_t ix=0; ix<xhits.GetEntries(); ix++)
221  {
222   AliSignal* sx=(AliSignal*)xhits.At(ix);
223   if (!sx) continue;
224   adcr=sx->GetSignal("ADC");
225   sigcor=fPe;
226   // If stored ADC data is un-calibrated, convert fPe to raw ADC value 
227   AliSignal* parent=(AliSignal*)sx->GetDevice();
228   if (parent)
229   {
230    if (parent->GetCalFunction("ADC"));
231    {
232     idrec=parent->GetUniqueID();
233     omr=(IceAOM*)fOmdb->GetObject(idrec,1);
234     TF1* fdecal=0;
235     if (omr) fdecal=omr->GetDecalFunction("ADC");
236     if (fdecal) sigcor=fdecal->Eval(fPe);
237    }
238   }
239   adcr=adcr-sigcor;
240   sx->SetSignal(adcr,"ADC");
241  }
242 }
243 ///////////////////////////////////////////////////////////////////////////