Josh, Irakli, David: AliCaloCalibSignal update - replacing possible 10k+ TProfiles...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRawUtils.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line 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. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17 /* History of cvs commits:
18  *
19  * $Log$
20  * Revision 1.10  2007/12/06 13:58:11  hristov
21  * Additional pritection. Do not delete the mapping, it is owned by another class
22  *
23  * Revision 1.9  2007/12/06 02:19:51  jklay
24  * incorporated fitting procedure from testbeam analysis into AliRoot
25  *
26  * Revision 1.8  2007/12/05 02:30:51  jklay
27  * modification to read Altro mappings into AliEMCALRecParam and pass to AliEMCALRawUtils from AliEMCALReconstructor; add option to AliEMCALRawUtils to set old RCU format (for testbeam) or not
28  *
29  * Revision 1.7  2007/11/14 15:51:46  gustavo
30  * Take out few unnecessary prints
31  *
32  * Revision 1.6  2007/11/01 01:23:51  mvl
33  * Removed call to SetOldRCUFormat, which is only needed for testbeam data
34  *
35  * Revision 1.5  2007/11/01 01:20:33  mvl
36  * Further improvement of peak finding; more robust fit
37  *
38  * Revision 1.4  2007/10/31 17:15:24  mvl
39  * Fixed bug in raw data unpacking; Added pedestal to signal fit; Added logic to deal with high/low gain
40  *
41  * Revision 1.3  2007/09/27 08:36:46  mvl
42  * More robust setting of fit range in FitRawSignal (P. Hristov)
43  *
44  * Revision 1.2  2007/09/03 20:55:35  jklay
45  * EMCAL e-by-e reconstruction methods from Cvetan
46  *
47  * Revision 1.1  2007/03/17 19:56:38  mvl
48  * Moved signal shape routines from AliEMCAL to separate class AliEMCALRawUtils to streamline raw data reconstruction code.
49  * */
50
51 //*-- Author: Marco van Leeuwen (LBL)
52 #include "AliEMCALRawUtils.h"
53
54 #include "TF1.h"
55 #include "TGraph.h"
56 #include "TSystem.h"
57
58 #include "AliLog.h"
59 #include "AliRun.h"
60 #include "AliRunLoader.h"
61 #include "AliCaloAltroMapping.h"
62 #include "AliAltroBuffer.h"
63 #include "AliRawReader.h"
64 #include "AliCaloRawStream.h"
65 #include "AliDAQ.h"
66
67 #include "AliEMCALRecParam.h"
68 #include "AliEMCALLoader.h"
69 #include "AliEMCALGeometry.h"
70 #include "AliEMCALDigitizer.h"
71 #include "AliEMCALDigit.h"
72 #include "AliEMCAL.h"
73
74 ClassImp(AliEMCALRawUtils)
75
76 // Signal shape parameters
77 Double_t AliEMCALRawUtils::fgTimeBinWidth  = 100E-9 ; // each sample is 100 ns
78 Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ;   // 15 time bins ~ 1.5 musec
79
80 // some digitization constants
81 Int_t    AliEMCALRawUtils::fgThreshold = 1;
82 Int_t    AliEMCALRawUtils::fgDDLPerSuperModule = 2;  // 2 ddls per SuperModule
83 Int_t    AliEMCALRawUtils::fgPedestalValue = 32;      // pedestal value for digits2raw
84 Double_t AliEMCALRawUtils::fgFEENoise = 3.;            // 3 ADC channels of noise (sampled)
85
86 AliEMCALRawUtils::AliEMCALRawUtils()
87   : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
88     fNPedSamples(0), fGeom(0), fOption("")
89 {
90
91   //These are default parameters.  
92   //Can be re-set from without with setter functions
93   fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits) 
94   fOrder = 2;                         // order of gamma fn
95   fTau = 2.35;                        // in units of timebin, from CERN 2007 testbeam
96   fNoiseThreshold = 3;
97   fNPedSamples = 5;
98
99   //Get Mapping RCU files from the AliEMCALRecParam                                 
100   const TObjArray* maps = AliEMCALRecParam::GetMappings();
101   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
102
103   for(Int_t i = 0; i < 2; i++) {
104     fMapping[i] = (AliAltroMapping*)maps->At(i);
105   }
106
107   //To make sure we match with the geometry in a simulation file,
108   //let's try to get it first.  If not, take the default geometry
109   AliRunLoader *rl = AliRunLoader::GetRunLoader();
110   if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) {
111     fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
112   } else {
113     AliInfo(Form("Using default geometry in raw reco"));
114     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
115   }
116
117   if(!fGeom) AliFatal(Form("Could not get geometry!"));
118
119 }
120
121 //____________________________________________________________________________
122 AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry)
123   : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
124     fNPedSamples(0), fGeom(pGeometry), fOption("")
125 {
126   //
127   // Initialize with the given geometry - constructor required by HLT
128   // HLT does not use/support AliRunLoader(s) instances
129   // This is a minimum intervention solution
130   // Comment by MPloskon@lbl.gov
131   //
132
133   //These are default parameters. 
134   //Can be re-set from without with setter functions 
135   fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits)
136   fOrder = 2;                         // order of gamma fn
137   fTau = 2.35;                        // in units of timebin, from CERN 2007 testbeam
138   fNoiseThreshold = 3;
139   fNPedSamples = 5;
140
141   //Get Mapping RCU files from the AliEMCALRecParam
142   const TObjArray* maps = AliEMCALRecParam::GetMappings();
143   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
144
145   for(Int_t i = 0; i < 2; i++) {
146     fMapping[i] = (AliAltroMapping*)maps->At(i);
147   }
148
149   if(!fGeom) AliFatal(Form("Could not get geometry!"));
150
151 }
152
153 //____________________________________________________________________________
154 AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU)
155   : TObject(),
156     fHighLowGainFactor(rawU.fHighLowGainFactor), 
157     fOrder(rawU.fOrder),
158     fTau(rawU.fTau),
159     fNoiseThreshold(rawU.fNoiseThreshold),
160     fNPedSamples(rawU.fNPedSamples),
161     fGeom(rawU.fGeom), 
162     fOption(rawU.fOption)
163 {
164   //copy ctor
165   fMapping[0] = rawU.fMapping[0];
166   fMapping[1] = rawU.fMapping[1];
167 }
168
169 //____________________________________________________________________________
170 AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU)
171 {
172   //assignment operator
173
174   if(this != &rawU) {
175     fHighLowGainFactor = rawU.fHighLowGainFactor;
176     fOrder = rawU.fOrder;
177     fTau = rawU.fTau;
178     fNoiseThreshold = rawU.fNoiseThreshold;
179     fNPedSamples = rawU.fNPedSamples;
180     fGeom = rawU.fGeom;
181     fOption = rawU.fOption;
182     fMapping[0] = rawU.fMapping[0];
183     fMapping[1] = rawU.fMapping[1];
184   }
185
186   return *this;
187
188 }
189
190 //____________________________________________________________________________
191 AliEMCALRawUtils::~AliEMCALRawUtils() {
192
193 }
194
195 //____________________________________________________________________________
196 void AliEMCALRawUtils::Digits2Raw()
197 {
198   // convert digits of the current event to raw data
199   
200   AliRunLoader *rl = AliRunLoader::GetRunLoader();
201   AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
202
203   // get the digits
204   loader->LoadDigits("EMCAL");
205   loader->GetEvent();
206   TClonesArray* digits = loader->Digits() ;
207   
208   if (!digits) {
209     Warning("Digits2Raw", "no digits found !");
210     return;
211   }
212
213   static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
214   AliAltroBuffer* buffers[nDDL];
215   for (Int_t i=0; i < nDDL; i++)
216     buffers[i] = 0;
217
218   Int_t adcValuesLow[fgkTimeBins];
219   Int_t adcValuesHigh[fgkTimeBins];
220
221   // loop over digits (assume ordered digits)
222   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
223     AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
224     if (digit->GetAmp() < fgThreshold) 
225       continue;
226
227     //get cell indices
228     Int_t nSM = 0;
229     Int_t nIphi = 0;
230     Int_t nIeta = 0;
231     Int_t iphi = 0;
232     Int_t ieta = 0;
233     Int_t nModule = 0;
234     fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
235     fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
236     
237     //Check which is the RCU of the cell.
238     Int_t iRCU = -111;
239     //RCU0
240     if (0<=iphi&&iphi<8) iRCU=0; // first cable row
241     else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; 
242     //second cable row
243     //RCU1
244     else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; 
245     //second cable row
246     else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
247     if (iRCU<0) 
248       Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU);
249     
250     //Which DDL?
251     Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU;
252     if (iDDL >= nDDL)
253       Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
254
255     if (buffers[iDDL] == 0) {      
256       // open new file and write dummy header
257       TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
258       buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCU]);
259       buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
260     }
261     
262     // out of time range signal (?)
263     if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
264       AliInfo("Signal is out of time range.\n");
265       buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp());
266       buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() );  // time bin
267       buffers[iDDL]->FillBuffer(3);          // bunch length      
268       buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM);  // trailer
269       // calculate the time response function
270     } else {
271       Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; 
272       if (lowgain) 
273         buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold);
274       else 
275         buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold);
276     }
277   }
278   
279   // write headers and close files
280   for (Int_t i=0; i < nDDL; i++) {
281     if (buffers[i]) {
282       buffers[i]->Flush();
283       buffers[i]->WriteDataHeader(kFALSE, kFALSE);
284       delete buffers[i];
285     }
286   }
287
288   loader->UnloadDigits();
289 }
290
291 //____________________________________________________________________________
292 void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
293 {
294   // convert raw data of the current event to digits                                                                                     
295
296   digitsArr->Clear(); 
297
298   if (!digitsArr) {
299     Error("Raw2Digits", "no digits found !");
300     return;
301   }
302   if (!reader) {
303     Error("Raw2Digits", "no raw reader found !");
304     return;
305   }
306
307   AliCaloRawStream in(reader,"EMCAL",fMapping);
308   // Select EMCAL DDL's;
309   reader->Select("EMCAL");
310
311   //Updated fitting routine from 2007 beam test takes into account
312   //possibility of two peaks in data and selects first one for fitting
313   //Also sets some of the starting parameters based on the shape of the
314   //given raw signal being fit
315
316   TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
317   signalF->SetParameters(10.,0.,fTau,fOrder,5.); //set all defaults once, just to be safe
318   signalF->SetParNames("amp","t0","tau","N","ped");
319   signalF->SetParameter(2,fTau); // tau in units of time bin
320   signalF->SetParLimits(2,2,-1);
321   signalF->SetParameter(3,fOrder); // order
322   signalF->SetParLimits(3,2,-1);
323   
324   Int_t id =  -1;
325   Float_t time = 0. ; 
326   Float_t amp = 0. ; 
327
328   //Graph to hold data we will fit (should be converted to an array
329   //later to speed up processing
330   TGraph * gSig = new TGraph(GetRawFormatTimeBins()); 
331
332   Int_t readOk = 1;
333   Int_t lowGain = 0;
334
335   while (readOk && in.GetModule() < 0) 
336     readOk = in.Next();  // Go to first digit
337
338   Int_t col = 0;
339   Int_t row = 0;
340
341   while (readOk) { 
342     id =  fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
343     lowGain = in.IsLowGain();
344     Int_t maxTime = in.GetTime();  // timebins come in reverse order
345     if (maxTime < 0 || maxTime >= GetRawFormatTimeBins()) {
346       AliWarning(Form("Invalid time bin %d",maxTime));
347       maxTime = GetRawFormatTimeBins();
348     }
349     gSig->Set(maxTime+1);
350     // There is some kind of zero-suppression in the raw data, 
351     // so set up the TGraph in advance
352     for (Int_t i=0; i < maxTime; i++) {
353       gSig->SetPoint(i, i , 0);
354     }
355
356     Int_t iTime = 0;
357     do {
358       if (in.GetTime() >= gSig->GetN()) {
359           AliWarning("Too many time bins");
360           gSig->Set(in.GetTime());
361       }
362       col = in.GetColumn();
363       row = in.GetRow();
364       
365       gSig->SetPoint(in.GetTime(), in.GetTime(), in.GetSignal()) ;
366
367       if (in.GetTime() > maxTime)
368         maxTime = in.GetTime();
369       iTime++;
370     } while ((readOk = in.Next()) && !in.IsNewHWAddress());
371
372     FitRaw(gSig, signalF, amp, time) ; 
373     
374     if (amp > 0 && amp < 10000) {  //check both high and low end of
375                                    //result, 10000 is somewhat arbitrary
376       AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
377       //cout << "col " << col-40 << " row " << row-8 << " lowGain " << lowGain << " amp " << amp << endl;
378
379       AddDigit(digitsArr, id, lowGain, (Int_t)amp, time);
380     }
381         
382     // Reset graph
383     for (Int_t index = 0; index < gSig->GetN(); index++) {
384       gSig->SetPoint(index, index, 0) ;  
385     } 
386     // Reset starting parameters for fit function
387     signalF->SetParameters(10.,0.,fTau,fOrder,5.); //reset all defaults just to be safe
388
389   }; // EMCAL entries loop
390   
391   delete signalF ; 
392   delete gSig;
393   
394   return ; 
395 }
396
397 //____________________________________________________________________________ 
398 void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) {
399   //
400   // Add a new digit. 
401   // This routine checks whether a digit exists already for this tower 
402   // and then decides whether to use the high or low gain info
403   //
404   // Called by Raw2Digits
405   
406   AliEMCALDigit *digit = 0, *tmpdigit = 0;
407   
408   TIter nextdigit(digitsArr);
409   while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) {
410     if (tmpdigit->GetId() == id)
411       digit = tmpdigit;
412   }
413
414   if (!digit) { // no digit existed for this tower; create one
415     if (lowGain) 
416       amp = Int_t(fHighLowGainFactor * amp); 
417     Int_t idigit = digitsArr->GetEntries();
418     new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ;   
419   }
420   else { // a digit already exists, check range 
421          // (use high gain if signal < cut value, otherwise low gain)
422     if (lowGain) { // new digit is low gain
423       if (digit->GetAmp() > fgkOverflowCut) {  // use if stored digit is out of range
424         digit->SetAmp(Int_t(fHighLowGainFactor * amp));
425         digit->SetTime(time);
426       }
427     }
428     else if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range
429       digit->SetAmp(amp);
430       digit->SetTime(time);
431     }
432   }
433 }
434
435 //____________________________________________________________________________ 
436 void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_t & time)
437 {
438   // Fits the raw signal time distribution; from AliEMCALGetter 
439
440   amp = time = 0. ; 
441   Double_t ped = 0;
442   Int_t nPed = 0;
443
444   for (Int_t index = 0; index < fNPedSamples; index++) {
445     Double_t ttime, signal;
446     gSig->GetPoint(index, ttime, signal) ; 
447     if (signal > 0) {
448       ped += signal;
449       nPed++;
450     }
451   }
452
453   if (nPed > 0)
454     ped /= nPed;
455   else {
456     AliWarning("Could not determine pedestal");   
457     ped = 10; // put some small value as first guess
458   }
459
460   Int_t max_found = 0;
461   Int_t i_max = 0;
462   Float_t max = -1;
463   Float_t max_fit = gSig->GetN();
464   Float_t min_after_sig = 9999;
465   Int_t tmin_after_sig = gSig->GetN();
466   Int_t n_ped_after_sig = 0;
467   Int_t plateau_width = 0;
468   Int_t plateau_start = 9999;
469
470   for (Int_t i=fNPedSamples; i < gSig->GetN(); i++) {
471     Double_t ttime, signal;
472     gSig->GetPoint(i, ttime, signal) ; 
473     if (!max_found && signal > max) {
474       i_max = i;
475       max = signal;
476     }
477     else if ( max > ped + fNoiseThreshold ) {
478       max_found = 1;
479       min_after_sig = signal;
480       tmin_after_sig = i;
481     }
482     if (max_found) {
483       if ( signal < min_after_sig) {
484         min_after_sig = signal;
485         tmin_after_sig = i;
486       }
487       if (i > tmin_after_sig + 5) {  // Two close peaks; end fit at minimum
488         max_fit = tmin_after_sig;
489         break;
490       }
491       if ( signal < ped + fNoiseThreshold)
492         n_ped_after_sig++;
493       if (n_ped_after_sig >= 5) {  // include 5 pedestal bins after peak
494         max_fit = i;
495         break;
496       }
497     }
498     //Add check on plateau
499     if (signal >= fgkRawSignalOverflow - fNoiseThreshold) {
500       if(plateau_width == 0) plateau_start = i;
501       plateau_width++;
502     }
503   }
504
505   if(plateau_width > 0) {
506     for(int j = 0; j < plateau_width-2; j++) {
507       //Note, have to remove the same point N times because after each
508       //remove, the positions of all subsequent points have shifted down
509       //We leave the first and last as anchor points
510       gSig->RemovePoint(plateau_start+1);
511     }
512   }
513
514   if ( max - ped > fNoiseThreshold ) { // else its noise 
515     AliDebug(2,Form("Fitting max %d ped %d", max, ped));
516     signalF->SetRange(0,max_fit);
517
518     if(max-ped > 50) 
519       signalF->SetParLimits(2,1,3);
520
521     signalF->SetParameter(4, ped) ; 
522     signalF->SetParameter(1, i_max);
523     signalF->SetParameter(0, max);
524     
525     gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
526     amp = signalF->GetParameter(0); 
527     time = signalF->GetParameter(1)*GetRawFormatTimeBinWidth() - fgTimeTrigger;
528   }
529   return;
530 }
531 //__________________________________________________________________
532 Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
533 {
534   // Matches version used in 2007 beam test
535   //
536   // Shape of the electronics raw reponse:
537   // It is a semi-gaussian, 2nd order Gamma function of the general form
538   //
539   // t' = (t - t0 + tau) / tau
540   // F = A * t**N * exp( N * ( 1 - t) )   for t >= 0
541   // F = 0                                for t < 0 
542   //
543   // parameters:
544   // A:   par[0]   // Amplitude = peak value
545   // t0:  par[1]
546   // tau: par[2]
547   // N:   par[3]
548   // ped: par[4]
549   //
550   Double_t signal ;
551   Double_t tau =par[2];
552   Double_t N =par[3];
553   Double_t ped = par[4];
554   Double_t xx = ( x[0] - par[1] + tau ) / tau ;
555
556   if (xx <= 0) 
557     signal = ped ;  
558   else {  
559     signal = ped + par[0] * TMath::Power(xx , N) * TMath::Exp(N * (1 - xx )) ; 
560   }
561   return signal ;  
562 }
563
564 //__________________________________________________________________
565 Bool_t AliEMCALRawUtils::RawSampledResponse(
566 const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const 
567 {
568   // for a start time dtime and an amplitude damp given by digit, 
569   // calculates the raw sampled response AliEMCAL::RawResponseFunction
570
571   Bool_t lowGain = kFALSE ; 
572
573   // A:   par[0]   // Amplitude = peak value
574   // t0:  par[1]                            
575   // tau: par[2]                            
576   // N:   par[3]                            
577   // ped: par[4]
578
579   TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 5);
580   signalF.SetParameter(0, damp) ; 
581   signalF.SetParameter(1, dtime + fgTimeTrigger) ; 
582   signalF.SetParameter(2, fTau) ; 
583   signalF.SetParameter(3, fOrder);
584   signalF.SetParameter(4, fgPedestalValue);
585
586   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
587     Double_t time = iTime * GetRawFormatTimeBinWidth() ;
588     Double_t signal = signalF.Eval(time) ;     
589
590     //According to Terry Awes, 13-Apr-2008
591     //add gaussian noise in quadrature to each sample
592     //Double_t noise = gRandom->Gaus(0.,fgFEENoise);
593     //signal = sqrt(signal*signal + noise*noise);
594
595     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
596     if ( adcH[iTime] > fgkRawSignalOverflow ){  // larger than 10 bits 
597       adcH[iTime] = fgkRawSignalOverflow ;
598       lowGain = kTRUE ; 
599     }
600
601     signal /= fHighLowGainFactor;
602
603     adcL[iTime] =  static_cast<Int_t>(signal + 0.5) ;
604     if ( adcL[iTime] > fgkRawSignalOverflow)  // larger than 10 bits 
605       adcL[iTime] = fgkRawSignalOverflow ;
606   }
607   return lowGain ; 
608 }