add SetReference method
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRawUtils.cxx
CommitLineData
ee299369 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$ */
ee299369 17
e5bbbc4e 18//_________________________________________________________________________
19// Utility Class for handling Raw data
20// Does all transitions from Digits to Raw and vice versa,
21// for simu and reconstruction
22//
23// Note: the current version is still simplified. Only
24// one raw signal per digit is generated; either high-gain or low-gain
25// Need to add concurrent high and low-gain info in the future
26// No pedestal is added to the raw signal.
ee299369 27//*-- Author: Marco van Leeuwen (LBL)
e5bbbc4e 28
ee299369 29#include "AliEMCALRawUtils.h"
507751ce 30#include <stdexcept>
21cad85c 31
4fe71e02 32#include "TF1.h"
33#include "TGraph.h"
6751f762 34#include <TRandom.h>
e5bbbc4e 35class TSystem;
21cad85c 36
e5bbbc4e 37class AliLog;
72c58de0 38#include "AliRun.h"
ee299369 39#include "AliRunLoader.h"
e5bbbc4e 40class AliCaloAltroMapping;
ee299369 41#include "AliAltroBuffer.h"
42#include "AliRawReader.h"
32cd4c24 43#include "AliCaloRawStreamV3.h"
ee299369 44#include "AliDAQ.h"
21cad85c 45
feedcab9 46#include "AliEMCALRecParam.h"
ee299369 47#include "AliEMCALLoader.h"
48#include "AliEMCALGeometry.h"
e5bbbc4e 49class AliEMCALDigitizer;
ee299369 50#include "AliEMCALDigit.h"
916f1e76 51#include "AliEMCALRawDigit.h"
20b636fc 52#include "AliEMCAL.h"
5e3106bc 53#include "AliCaloCalibPedestal.h"
9f467289 54#include "AliCaloFastAltroFitv0.h"
c8603a2b 55#include "AliCaloNeuralFit.h"
16605c06 56#include "AliCaloBunchInfo.h"
57#include "AliCaloFitResults.h"
7683df1d 58#include "AliCaloRawAnalyzerFastFit.h"
59#include "AliCaloRawAnalyzerNN.h"
16605c06 60#include "AliCaloRawAnalyzerLMS.h"
61#include "AliCaloRawAnalyzerPeakFinder.h"
62#include "AliCaloRawAnalyzerCrude.h"
9f467289 63
ee299369 64ClassImp(AliEMCALRawUtils)
21cad85c 65
ee299369 66// Signal shape parameters
89d338a6 67Int_t AliEMCALRawUtils::fgTimeBins = 256; // number of sampling bins of the raw RO signal (we typically use 15-50; theoretical max is 1k+)
e5bbbc4e 68Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns
09974781 69Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ; // 15 time bins ~ 1.5 musec
ee299369 70
71// some digitization constants
72Int_t AliEMCALRawUtils::fgThreshold = 1;
73Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule
916f1e76 74Int_t AliEMCALRawUtils::fgPedestalValue = 0; // pedestal value for digits2raw, default generate ZS data
e5bbbc4e 75Double_t AliEMCALRawUtils::fgFEENoise = 3.; // 3 ADC channels of noise (sampled)
ee299369 76
16605c06 77AliEMCALRawUtils::AliEMCALRawUtils(fitAlgorithm fitAlgo)
b4133f05 78 : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
9f467289 79 fNPedSamples(0), fGeom(0), fOption(""),
f4a4dc82 80 fRemoveBadChannels(kTRUE),fFittingAlgorithm(0),
81 fTimeMin(-1.),fTimeMax(1.),
82 fUseFALTRO(kFALSE),fRawAnalyzer(0)
8cb998bd 83{
b4133f05 84
85 //These are default parameters.
86 //Can be re-set from without with setter functions
9f467289 87 //Already set in the OCDB and passed via setter in the AliEMCALReconstructor
46f1d25f 88 fHighLowGainFactor = 16. ; // Adjusted for a low gain range of 82 GeV (10 bits)
89 fOrder = 2; // Order of gamma fn
90 fTau = 2.35; // in units of timebin, from CERN 2007 testbeam
91 fNoiseThreshold = 3; // 3 ADC counts is approx. noise level
92 fNPedSamples = 4; // Less than this value => likely pedestal samples
93 fRemoveBadChannels = kFALSE; // Do not remove bad channels before fitting
94 fUseFALTRO = kTRUE; // Get the trigger FALTRO information and pass it to digits.
4fe71e02 95 SetFittingAlgorithm(fitAlgo);
16605c06 96
65bdc82f 97 //Get Mapping RCU files from the AliEMCALRecParam
98 const TObjArray* maps = AliEMCALRecParam::GetMappings();
99 if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
100
21cad85c 101 for(Int_t i = 0; i < 4; i++) {
65bdc82f 102 fMapping[i] = (AliAltroMapping*)maps->At(i);
103 }
104
72c58de0 105 //To make sure we match with the geometry in a simulation file,
106 //let's try to get it first. If not, take the default geometry
33c3c91a 107 AliRunLoader *rl = AliRunLoader::Instance();
916f1e76 108 if (rl && rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) {
72c58de0 109 fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
110 } else {
111 AliInfo(Form("Using default geometry in raw reco"));
937d0661 112 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
65bdc82f 113 }
114
72c58de0 115 if(!fGeom) AliFatal(Form("Could not get geometry!"));
116
65bdc82f 117}
118
119//____________________________________________________________________________
16605c06 120AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry, fitAlgorithm fitAlgo)
5544799a 121 : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
9f467289 122 fNPedSamples(0), fGeom(pGeometry), fOption(""),
f4a4dc82 123 fRemoveBadChannels(kTRUE),fFittingAlgorithm(0),
124 fTimeMin(-1.),fTimeMax(1.),
125 fUseFALTRO(kFALSE),fRawAnalyzer()
5544799a 126{
127 //
128 // Initialize with the given geometry - constructor required by HLT
129 // HLT does not use/support AliRunLoader(s) instances
130 // This is a minimum intervention solution
131 // Comment by MPloskon@lbl.gov
132 //
133
134 //These are default parameters.
135 //Can be re-set from without with setter functions
9f467289 136 //Already set in the OCDB and passed via setter in the AliEMCALReconstructor
46f1d25f 137 fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits)
138 fOrder = 2; // order of gamma fn
139 fTau = 2.35; // in units of timebin, from CERN 2007 testbeam
140 fNoiseThreshold = 3; // 3 ADC counts is approx. noise level
141 fNPedSamples = 4; // Less than this value => likely pedestal samples
142 fRemoveBadChannels = kFALSE; // Do not remove bad channels before fitting
143 fUseFALTRO = kTRUE; // Get the trigger FALTRO information and pass it to digits.
4fe71e02 144 SetFittingAlgorithm(fitAlgo);
46f1d25f 145
5544799a 146 //Get Mapping RCU files from the AliEMCALRecParam
147 const TObjArray* maps = AliEMCALRecParam::GetMappings();
148 if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
149
21cad85c 150 for(Int_t i = 0; i < 4; i++) {
5544799a 151 fMapping[i] = (AliAltroMapping*)maps->At(i);
152 }
153
154 if(!fGeom) AliFatal(Form("Could not get geometry!"));
155
156}
157
158//____________________________________________________________________________
65bdc82f 159AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU)
160 : TObject(),
161 fHighLowGainFactor(rawU.fHighLowGainFactor),
b4133f05 162 fOrder(rawU.fOrder),
163 fTau(rawU.fTau),
164 fNoiseThreshold(rawU.fNoiseThreshold),
165 fNPedSamples(rawU.fNPedSamples),
65bdc82f 166 fGeom(rawU.fGeom),
9f467289 167 fOption(rawU.fOption),
168 fRemoveBadChannels(rawU.fRemoveBadChannels),
16605c06 169 fFittingAlgorithm(rawU.fFittingAlgorithm),
f4a4dc82 170 fTimeMin(rawU.fTimeMin),fTimeMax(rawU.fTimeMax),
46f1d25f 171 fUseFALTRO(rawU.fUseFALTRO),
16605c06 172 fRawAnalyzer(rawU.fRawAnalyzer)
65bdc82f 173{
174 //copy ctor
175 fMapping[0] = rawU.fMapping[0];
176 fMapping[1] = rawU.fMapping[1];
21cad85c 177 fMapping[2] = rawU.fMapping[2];
178 fMapping[3] = rawU.fMapping[3];
65bdc82f 179}
180
181//____________________________________________________________________________
182AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU)
183{
184 //assignment operator
185
186 if(this != &rawU) {
187 fHighLowGainFactor = rawU.fHighLowGainFactor;
46f1d25f 188 fOrder = rawU.fOrder;
189 fTau = rawU.fTau;
190 fNoiseThreshold = rawU.fNoiseThreshold;
191 fNPedSamples = rawU.fNPedSamples;
192 fGeom = rawU.fGeom;
193 fOption = rawU.fOption;
9f467289 194 fRemoveBadChannels = rawU.fRemoveBadChannels;
195 fFittingAlgorithm = rawU.fFittingAlgorithm;
f4a4dc82 196 fTimeMin = rawU.fTimeMin;
197 fTimeMax = rawU.fTimeMax;
46f1d25f 198 fUseFALTRO = rawU.fUseFALTRO;
199 fRawAnalyzer = rawU.fRawAnalyzer;
200 fMapping[0] = rawU.fMapping[0];
201 fMapping[1] = rawU.fMapping[1];
202 fMapping[2] = rawU.fMapping[2];
203 fMapping[3] = rawU.fMapping[3];
65bdc82f 204 }
205
206 return *this;
207
ee299369 208}
65bdc82f 209
ee299369 210//____________________________________________________________________________
211AliEMCALRawUtils::~AliEMCALRawUtils() {
e5bbbc4e 212 //dtor
65bdc82f 213
ee299369 214}
65bdc82f 215
ee299369 216//____________________________________________________________________________
65bdc82f 217void AliEMCALRawUtils::Digits2Raw()
ee299369 218{
219 // convert digits of the current event to raw data
220
33c3c91a 221 AliRunLoader *rl = AliRunLoader::Instance();
ee299369 222 AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
223
224 // get the digits
225 loader->LoadDigits("EMCAL");
226 loader->GetEvent();
227 TClonesArray* digits = loader->Digits() ;
228
229 if (!digits) {
230 Warning("Digits2Raw", "no digits found !");
231 return;
232 }
65bdc82f 233
ee299369 234 static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
235 AliAltroBuffer* buffers[nDDL];
236 for (Int_t i=0; i < nDDL; i++)
237 buffers[i] = 0;
238
e2c2134b 239 TArrayI adcValuesLow(fgTimeBins);
240 TArrayI adcValuesHigh(fgTimeBins);
ee299369 241
ee299369 242 // loop over digits (assume ordered digits)
243 for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
244 AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
829ba234 245 if (digit->GetAmplitude() < fgThreshold)
ee299369 246 continue;
247
248 //get cell indices
249 Int_t nSM = 0;
250 Int_t nIphi = 0;
251 Int_t nIeta = 0;
252 Int_t iphi = 0;
253 Int_t ieta = 0;
254 Int_t nModule = 0;
65bdc82f 255 fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
256 fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
ee299369 257
21cad85c 258 //Check which is the RCU, 0 or 1, of the cell.
ee299369 259 Int_t iRCU = -111;
260 //RCU0
261 if (0<=iphi&&iphi<8) iRCU=0; // first cable row
262 else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half;
263 //second cable row
264 //RCU1
265 else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half;
266 //second cable row
267 else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
21cad85c 268
269 if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
270
e36e3bcf 271 if (iRCU<0)
272 Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU);
ee299369 273
274 //Which DDL?
275 Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU;
276 if (iDDL >= nDDL)
277 Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
278
279 if (buffers[iDDL] == 0) {
280 // open new file and write dummy header
281 TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
21cad85c 282 //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C
283 Int_t iRCUside=iRCU+(nSM%2)*2;
284 //iRCU=0 and even (0) SM -> RCU0A.data 0
285 //iRCU=1 and even (0) SM -> RCU1A.data 1
286 //iRCU=0 and odd (1) SM -> RCU0C.data 2
287 //iRCU=1 and odd (1) SM -> RCU1C.data 3
288 //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl;
289 buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]);
ee299369 290 buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy;
291 }
292
293 // out of time range signal (?)
294 if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
295 AliInfo("Signal is out of time range.\n");
829ba234 296 buffers[iDDL]->FillBuffer((Int_t)digit->GetAmplitude());
ee299369 297 buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin
298 buffers[iDDL]->FillBuffer(3); // bunch length
299 buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer
300 // calculate the time response function
301 } else {
829ba234 302 Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ;
ee299369 303 if (lowgain)
e2c2134b 304 buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow.GetArray(), fgThreshold);
ee299369 305 else
e2c2134b 306 buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh.GetArray(), fgThreshold);
ee299369 307 }
308 }
309
310 // write headers and close files
311 for (Int_t i=0; i < nDDL; i++) {
312 if (buffers[i]) {
313 buffers[i]->Flush();
314 buffers[i]->WriteDataHeader(kFALSE, kFALSE);
315 delete buffers[i];
316 }
317 }
65bdc82f 318
ee299369 319 loader->UnloadDigits();
320}
321
322//____________________________________________________________________________
916f1e76 323void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap, TClonesArray *digitsTRG)
ee299369 324{
65bdc82f 325 // convert raw data of the current event to digits
ee299369 326
c47157cd 327 digitsArr->Clear();
ee299369 328
c47157cd 329 if (!digitsArr) {
ee299369 330 Error("Raw2Digits", "no digits found !");
331 return;
332 }
333 if (!reader) {
334 Error("Raw2Digits", "no raw reader found !");
335 return;
336 }
337
32cd4c24 338 AliCaloRawStreamV3 in(reader,"EMCAL",fMapping);
ee299369 339 // Select EMCAL DDL's;
7643e728 340 reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL
feedcab9 341
16605c06 342 // fRawAnalyzer setup
f57baa2d 343 fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done
16605c06 344 fRawAnalyzer->SetAmpCut(fNoiseThreshold);
345 fRawAnalyzer->SetFitArrayCut(fNoiseThreshold);
346 fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later
ee299369 347
16605c06 348 // channel info parameters
829ba234 349 Int_t lowGain = 0;
e5bbbc4e 350 Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref.
ee299369 351
32cd4c24 352 // start loop over input stream
353 while (in.NextDDL()) {
916f1e76 354
355// if ( in.GetDDLNumber() != 0 && in.GetDDLNumber() != 2 ) continue;
356
32cd4c24 357 while (in.NextChannel()) {
7643e728 358
916f1e76 359/*
360 Int_t hhwAdd = in.GetHWAddress();
361 UShort_t iiBranch = ( hhwAdd >> 11 ) & 0x1; // 0/1
362 UShort_t iiFEC = ( hhwAdd >> 7 ) & 0xF;
363 UShort_t iiChip = ( hhwAdd >> 4 ) & 0x7;
364 UShort_t iiChannel = hhwAdd & 0xF;
365
366 if ( !( iiBranch == 0 && iiFEC == 1 && iiChip == 3 && ( iiChannel >= 8 && iiChannel <= 15 ) ) && !( iiBranch == 1 && iiFEC == 0 && in.GetColumn() == 0 ) ) continue;
367*/
368
7643e728 369 //Check if the signal is high or low gain and then do the fit,
16605c06 370 //if it is from TRU or LEDMon do not fit
7643e728 371 caloFlag = in.GetCaloFlag();
916f1e76 372// if (caloFlag != 0 && caloFlag != 1) continue;
373 if (caloFlag > 2) continue; // Work with ALTRO and FALTRO
374
375 //Do not fit bad channels of ALTRO
376 if(caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) {
5e3106bc 377 //printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow());
378 continue;
379 }
380
16605c06 381 vector<AliCaloBunchInfo> bunchlist;
32cd4c24 382 while (in.NextBunch()) {
16605c06 383 bunchlist.push_back( AliCaloBunchInfo(in.GetStartTimeBin(), in.GetBunchLength(), in.GetSignals() ) );
384 } // loop over bunches
7643e728 385
916f1e76 386
507751ce 387 if ( caloFlag < 2 ){ // ALTRO
916f1e76 388
507751ce 389 Float_t time = 0;
829ba234 390 Float_t amp = 0;
391 short timeEstimate = 0;
507751ce 392 Float_t ampEstimate = 0;
393 Bool_t fitDone = kFALSE;
916f1e76 394
7683df1d 395 if ( fFittingAlgorithm == kFastFit || fFittingAlgorithm == kNeuralNet || fFittingAlgorithm == kLMS || fFittingAlgorithm == kPeakFinder || fFittingAlgorithm == kCrude) {
16605c06 396 // all functionality to determine amp and time etc is encapsulated inside the Evaluate call for these methods
397 AliCaloFitResults fitResults = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2());
398
829ba234 399 amp = fitResults.GetAmp();
400 time = fitResults.GetTime();
507751ce 401 timeEstimate = fitResults.GetMaxTimebin();
829ba234 402 ampEstimate = fitResults.GetMaxSig();
507751ce 403 if (fitResults.GetStatus() == AliCaloFitResults::kFitPar) {
404 fitDone = kTRUE;
405 }
16605c06 406 }
407 else { // for the other methods we for now use the functionality of
408 // AliCaloRawAnalyzer as well, to select samples and prepare for fits,
409 // if it looks like there is something to fit
410
411 // parameters init.
507751ce 412 Float_t pedEstimate = 0;
16605c06 413 short maxADC = 0;
16605c06 414 Int_t first = 0;
415 Int_t last = 0;
416 Int_t bunchIndex = 0;
417 //
418 // The PreFitEvaluateSamples + later call to FitRaw will hopefully
419 // be replaced by a single Evaluate call or so soon, like for the other
420 // methods, but this should be good enough for evaluation of
421 // the methods for now (Jan. 2010)
422 //
423 int nsamples = fRawAnalyzer->PreFitEvaluateSamples( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2(), bunchIndex, ampEstimate, maxADC, timeEstimate, pedEstimate, first, last);
7643e728 424
16605c06 425 if (ampEstimate > fNoiseThreshold) { // something worth looking at
7643e728 426
f57baa2d 427 time = timeEstimate; // maxrev in AliCaloRawAnalyzer speak; comes with an offset w.r.t. real timebin
428 Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
16605c06 429 amp = ampEstimate;
430
431 if ( nsamples > 1 ) { // possibly something to fit
507751ce 432 FitRaw(first, last, amp, time, fitDone);
f57baa2d 433 time += timebinOffset;
507751ce 434 timeEstimate += timebinOffset;
9f467289 435 }
16605c06 436
16605c06 437 } // ampEstimate check
438 } // method selection
507751ce 439
440 if ( fitDone ) { // brief sanity check of fit results
441 Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
442 Float_t timeDiff = time - timeEstimate;
443 if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) ) {
444 // AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations amp %f time %d", amp, time, ampEstimate, timeEstimate));
445
446 // for now just overwrite the fit results with the simple/initial estimate
829ba234 447 amp = ampEstimate;
448 time = timeEstimate;
507751ce 449 fitDone = kFALSE;
450 }
451 } // fitDone
16605c06 452
029fe7a2 453 if (amp > fNoiseThreshold && amp<fgkRawSignalOverflow) { // something to be stored
507751ce 454 if ( ! fitDone) { // smear ADC with +- 0.5 uniform (avoid discrete effects)
455 amp += (0.5 - gRandom->Rndm()); // Rndm generates a number in ]0,1]
456 }
457
829ba234 458 Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
459 lowGain = in.IsLowGain();
7643e728 460
16605c06 461 // go from time-bin units to physical time fgtimetrigger
462 time = time * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger?
e7acbc37 463 // subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0:
464 time -= in.GetL1Phase();
7643e728 465
466 AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
467 // printf("Added tower: SM %d, row %d, column %d, amp %3.2f\n",in.GetModule(), in.GetRow(), in.GetColumn(),amp);
829ba234 468 AddDigit(digitsArr, id, lowGain, amp, time);
7643e728 469 }
470
916f1e76 471 }//ALTRO
46f1d25f 472 else if(fUseFALTRO)
916f1e76 473 {// Fake ALTRO
474 // if (maxTimeBin && gSig->GetN() > maxTimeBin + 10) gSig->Set(maxTimeBin + 10); // set actual max size of TGraph
916f1e76 475 Int_t hwAdd = in.GetHWAddress();
476 UShort_t iRCU = in.GetDDLNumber() % 2; // 0/1
477 UShort_t iBranch = ( hwAdd >> 11 ) & 0x1; // 0/1
478
479 // Now find TRU number
480 Int_t itru = 3 * in.GetModule() + ( (iRCU << 1) | iBranch ) - 1;
481
482 AliDebug(1,Form("Found TRG digit in TRU: %2d ADC: %2d",itru,in.GetColumn()));
483
484 Int_t idtrg;
485
486 Bool_t isOK = fGeom->GetAbsFastORIndexFromTRU(itru, in.GetColumn(), idtrg);
487
488 Int_t timeSamples[256]; for (Int_t j=0;j<256;j++) timeSamples[j] = 0;
489 Int_t nSamples = 0;
490
491 for (std::vector<AliCaloBunchInfo>::iterator itVectorData = bunchlist.begin(); itVectorData != bunchlist.end(); itVectorData++)
492 {
493 AliCaloBunchInfo bunch = *(itVectorData);
494
495 const UShort_t* sig = bunch.GetData();
496 Int_t startBin = bunch.GetStartBin();
497
498 for (Int_t iS = 0; iS < bunch.GetLength(); iS++)
499 {
500 Int_t time = startBin--;
501 Int_t amp = sig[iS];
502
503 if ( amp ) timeSamples[nSamples++] = ( ( time << 12 ) & 0xFF000 ) | ( amp & 0xFFF );
504 }
505 }
506
507 if (nSamples && isOK) AddDigit(digitsTRG, idtrg, timeSamples, nSamples);
508 }//Fake ALTRO
32cd4c24 509 } // end while over channel
510 } //end while over DDL's, of input stream
16605c06 511
f4a4dc82 512 TrimDigits(digitsArr);
513
ee299369 514 return ;
515}
516
517//____________________________________________________________________________
916f1e76 518void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t timeSamples[], Int_t nSamples)
519{
f4a4dc82 520 //Add raw sample to raw digit
521 new((*digitsArr)[digitsArr->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
522
523 // Int_t idx = digitsArr->GetEntriesFast()-1;
524 // AliEMCALRawDigit* d = (AliEMCALRawDigit*)digitsArr->At(idx);
916f1e76 525}
526
527//____________________________________________________________________________
829ba234 528void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time) {
82cbdfca 529 //
530 // Add a new digit.
531 // This routine checks whether a digit exists already for this tower
532 // and then decides whether to use the high or low gain info
533 //
534 // Called by Raw2Digits
535
536 AliEMCALDigit *digit = 0, *tmpdigit = 0;
82cbdfca 537 TIter nextdigit(digitsArr);
538 while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) {
f4a4dc82 539 if (tmpdigit->GetId() == id) digit = tmpdigit;
82cbdfca 540 }
541
542 if (!digit) { // no digit existed for this tower; create one
f4a4dc82 543 Int_t type = AliEMCALDigit::kHG; // use enum in AliEMCALDigit
544 if (lowGain) {
545 amp *= fHighLowGainFactor;
546 type = AliEMCALDigit::kLGnoHG;
547 }
548 Int_t idigit = digitsArr->GetEntries();
549 new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit) ;
550 AliDebug(2,Form("Add digit Id %d for the first time, type %d", id, type));
551 }//digit added first time
82cbdfca 552 else { // a digit already exists, check range
f4a4dc82 553 // (use high gain if signal < cut value, otherwise low gain)
554 if (lowGain) { // new digit is low gain
555 if (digit->GetAmplitude() > fgkOverflowCut) { // use if previously stored (HG) digit is out of range
556 digit->SetAmplitude(fHighLowGainFactor * amp);
557 digit->SetTime(time);
558 digit->SetType(AliEMCALDigit::kLG);
559 AliDebug(2,Form("Add LG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType()));
560 }
561 }//new low gain digit
562 else { // new digit is high gain
563 if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range
564 digit->SetAmplitude(amp);
565 digit->SetTime(time);
566 digit->SetType(AliEMCALDigit::kHG);
567 AliDebug(2,Form("Add HG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType()));
568 }
569 else { // HG out of range, just change flag value to show that HG did exist
570 digit->SetType(AliEMCALDigit::kLG);
571 AliDebug(2,Form("Change LG digit to HG, ID %d, type %d", digit->GetId(), digit->GetType()));
572 }
573 }//new high gain digit
574 }//digit existed replace it
575
576}
577
578//____________________________________________________________________________
579void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr)
580{
581 // Remove digits with only low gain and large time
582
583 AliEMCALDigit *digit = 0;
584 Int_t n = 0;
585 Int_t nDigits = digitsArr->GetEntriesFast();
586 TIter nextdigit(digitsArr);
587 while ((digit = (AliEMCALDigit*) nextdigit())) {
588
589 //Check if only LG existed, remove if so
590 if (digit->GetType() == AliEMCALDigit::kLGnoHG) {
591 AliDebug(1,Form("Remove digit with id %d, LGnoHG",digit->GetId()));
592 digitsArr->Remove(digit);
82cbdfca 593 }
f4a4dc82 594 //Check if time if too large or too small, remove if so
595 else if(fTimeMin > digit->GetTime() || fTimeMax < digit->GetTime()) {
596 digitsArr->Remove(digit);
597 AliDebug(1,Form("Remove digit with id %d, Bad Time %e",digit->GetId(), digit->GetTime()));
82cbdfca 598 }
f4a4dc82 599 //Good digit, just reassign the index of the digit in case there was a previous removal
600 else {
601 digit->SetIndexInList(n);
602 n++;
603 }
604 }//while
605
606 digitsArr->Compress();
607 AliDebug(1,Form("N Digits before trimming : %d; after array compression %d",nDigits,digitsArr->GetEntriesFast()));
608
82cbdfca 609}
f4a4dc82 610
82cbdfca 611//____________________________________________________________________________
507751ce 612void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Bool_t & fitDone) const
16605c06 613{ // Fits the raw signal time distribution
614
615 //--------------------------------------------------
616 //Do the fit, different fitting algorithms available
617 //--------------------------------------------------
618 int nsamples = lastTimeBin - firstTimeBin + 1;
507751ce 619 fitDone = kFALSE;
ee299369 620
16605c06 621 switch(fFittingAlgorithm) {
622 case kStandard:
623 {
7683df1d 624 if (nsamples < 3) { return; } // nothing much to fit
16605c06 625 //printf("Standard fitter \n");
7683df1d 626
16605c06 627 // Create Graph to hold data we will fit
7683df1d 628 TGraph *gSig = new TGraph( nsamples);
629 for (int i=0; i<nsamples; i++) {
630 Int_t timebin = firstTimeBin + i;
f57baa2d 631 gSig->SetPoint(i, timebin, fRawAnalyzer->GetReversed(timebin));
7683df1d 632 }
633
16605c06 634 TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
635 signalF->SetParameters(10.,5.,fTau,fOrder,0.); //set all defaults once, just to be safe
636 signalF->SetParNames("amp","t0","tau","N","ped");
637 signalF->FixParameter(2,fTau); // tau in units of time bin
638 signalF->FixParameter(3,fOrder); // order
639 signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here
640 signalF->SetParameter(1, time);
641 signalF->SetParameter(0, amp);
507751ce 642 // set rather loose parameter limits
643 signalF->SetParLimits(0, 0.5*amp, 2*amp );
644 signalF->SetParLimits(1, time - 4, time + 4);
645
646 try {
647 gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
648 // assign fit results
829ba234 649 amp = signalF->GetParameter(0);
507751ce 650 time = signalF->GetParameter(1);
651
652 // cross-check with ParabolaFit to see if the results make sense
653 FitParabola(gSig, amp); // amp is possibly updated
654 fitDone = kTRUE;
655 }
656 catch (const std::exception & e) {
657 AliError( Form("TGraph Fit exception %s", e.what()) );
658 // stay with default amp and time in case of exception, i.e. no special action required
659 fitDone = kFALSE;
660 }
16605c06 661 delete signalF;
662
16605c06 663 //printf("Std : Amp %f, time %g\n",amp, time);
7683df1d 664 delete gSig; // delete TGraph
16605c06 665
666 break;
667 }//kStandard Fitter
668 //----------------------------
7683df1d 669 case kLogFit:
16605c06 670 {
7683df1d 671 if (nsamples < 3) { return; } // nothing much to fit
672 //printf("LogFit \n");
673
674 // Create Graph to hold data we will fit
675 TGraph *gSigLog = new TGraph( nsamples);
676 for (int i=0; i<nsamples; i++) {
677 Int_t timebin = firstTimeBin + i;
678 gSigLog->SetPoint(timebin, timebin, TMath::Log(fRawAnalyzer->GetReversed(timebin) ) );
7643e728 679 }
7683df1d 680
681 TF1 * signalFLog = new TF1("signalLog", RawResponseFunctionLog, 0, GetRawFormatTimeBins(), 5);
682 signalFLog->SetParameters(2.3, 5.,fTau,fOrder,0.); //set all defaults once, just to be safe
683 signalFLog->SetParNames("amplog","t0","tau","N","ped");
684 signalFLog->FixParameter(2,fTau); // tau in units of time bin
685 signalFLog->FixParameter(3,fOrder); // order
686 signalFLog->FixParameter(4, 0); // pedestal should be subtracted when we get here
687 signalFLog->SetParameter(1, time);
688 if (amp>=1) {
689 signalFLog->SetParameter(0, TMath::Log(amp));
16605c06 690 }
7683df1d 691
692 gSigLog->Fit(signalFLog, "QROW"); // Note option 'W': equal errors on all points
693
694 // assign fit results
695 Double_t amplog = signalFLog->GetParameter(0); //Not Amp, but Log of Amp
696 amp = TMath::Exp(amplog);
697 time = signalFLog->GetParameter(1);
507751ce 698 fitDone = kTRUE;
7683df1d 699
700 delete signalFLog;
701 //printf("LogFit: Amp %f, time %g\n",amp, time);
702 delete gSigLog;
16605c06 703 break;
7683df1d 704 } //kLogFit
705 //----------------------------
706
16605c06 707 //----------------------------
708 }//switch fitting algorithms
fb070798 709
16605c06 710 return;
711}
8cb998bd 712
16605c06 713//__________________________________________________________________
714void AliEMCALRawUtils::FitParabola(const TGraph *gSig, Float_t & amp) const
715{
716 //BEG YS alternative methods to calculate the amplitude
717 Double_t * ymx = gSig->GetX() ;
718 Double_t * ymy = gSig->GetY() ;
719 const Int_t kN = 3 ;
720 Double_t ymMaxX[kN] = {0., 0., 0.} ;
721 Double_t ymMaxY[kN] = {0., 0., 0.} ;
722 Double_t ymax = 0. ;
723 // find the maximum amplitude
724 Int_t ymiMax = 0 ;
725 for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
726 if (ymy[ymi] > ymMaxY[0] ) {
727 ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
728 ymMaxX[0] = ymx[ymi] ;
729 ymiMax = ymi ;
730 }
731 }
732 // find the maximum by fitting a parabola through the max and the two adjacent samples
733 if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
734 ymMaxY[1] = ymy[ymiMax+1] ;
735 ymMaxY[2] = ymy[ymiMax-1] ;
736 ymMaxX[1] = ymx[ymiMax+1] ;
737 ymMaxX[2] = ymx[ymiMax-1] ;
738 if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
739 //fit a parabola through the 3 points y= a+bx+x*x*x
740 Double_t sy = 0 ;
741 Double_t sx = 0 ;
742 Double_t sx2 = 0 ;
743 Double_t sx3 = 0 ;
744 Double_t sx4 = 0 ;
745 Double_t sxy = 0 ;
746 Double_t sx2y = 0 ;
747 for (Int_t i = 0; i < kN ; i++) {
748 sy += ymMaxY[i] ;
749 sx += ymMaxX[i] ;
750 sx2 += ymMaxX[i]*ymMaxX[i] ;
751 sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
752 sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
753 sxy += ymMaxX[i]*ymMaxY[i] ;
754 sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ;
755 }
756 Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2);
757 Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
758 Double_t c = cN / cD ;
759 Double_t b = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
760 Double_t a = (sy-b*sx-c*sx2)/kN ;
761 Double_t xmax = -b/(2*c) ;
762 ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
029fe7a2 763 amp = ymax;
16605c06 764 }
765 }
766
767 Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ;
768 if (diff > 0.1)
769 amp = ymMaxY[0] ;
770 //printf("Yves : Amp %f, time %g\n",amp, time);
771 //END YS
ee299369 772 return;
773}
16605c06 774
ee299369 775//__________________________________________________________________
776Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
777{
8cb998bd 778 // Matches version used in 2007 beam test
779 //
ee299369 780 // Shape of the electronics raw reponse:
781 // It is a semi-gaussian, 2nd order Gamma function of the general form
782 //
7643e728 783 // xx = (t - t0 + tau) / tau [xx is just a convenient help variable]
784 // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0
785 // F = 0 for xx < 0
ee299369 786 //
787 // parameters:
8cb998bd 788 // A: par[0] // Amplitude = peak value
789 // t0: par[1]
790 // tau: par[2]
791 // N: par[3]
792 // ped: par[4]
ee299369 793 //
794 Double_t signal ;
8cb998bd 795 Double_t tau =par[2];
e5bbbc4e 796 Double_t n =par[3];
8cb998bd 797 Double_t ped = par[4];
798 Double_t xx = ( x[0] - par[1] + tau ) / tau ;
ee299369 799
5a056daa 800 if (xx <= 0)
8cb998bd 801 signal = ped ;
ee299369 802 else {
e5bbbc4e 803 signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
ee299369 804 }
805 return signal ;
806}
807
808//__________________________________________________________________
7683df1d 809Double_t AliEMCALRawUtils::RawResponseFunctionLog(Double_t *x, Double_t *par)
810{
811 // Matches version used in 2007 beam test
812 //
813 // Shape of the electronics raw reponse:
814 // It is a semi-gaussian, 2nd order Gamma function of the general form
815 //
816 // xx = (t - t0 + tau) / tau [xx is just a convenient help variable]
817 // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0
818 // F = 0 for xx < 0
819 //
820 // parameters:
821 // Log[A]: par[0] // Amplitude = peak value
822 // t0: par[1]
823 // tau: par[2]
824 // N: par[3]
825 // ped: par[4]
826 //
827 Double_t signal ;
828 Double_t tau =par[2];
829 Double_t n =par[3];
830 //Double_t ped = par[4]; // not used
831 Double_t xx = ( x[0] - par[1] + tau ) / tau ;
832
833 if (xx < 0)
834 signal = par[0] - n*TMath::Log(TMath::Abs(xx)) + n * (1 - xx ) ;
835 else {
836 signal = par[0] + n*TMath::Log(xx) + n * (1 - xx ) ;
837 }
838 return signal ;
839}
840
841//__________________________________________________________________
6751f762 842Bool_t AliEMCALRawUtils::RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL, const Int_t keyErr) const
ee299369 843{
844 // for a start time dtime and an amplitude damp given by digit,
845 // calculates the raw sampled response AliEMCAL::RawResponseFunction
846
ee299369 847 Bool_t lowGain = kFALSE ;
848
48a56166 849 // A: par[0] // Amplitude = peak value
850 // t0: par[1]
851 // tau: par[2]
852 // N: par[3]
853 // ped: par[4]
854
56e13066 855 TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
48a56166 856 signalF.SetParameter(0, damp) ;
56e13066 857 signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ;
b4133f05 858 signalF.SetParameter(2, fTau) ;
859 signalF.SetParameter(3, fOrder);
fe93d365 860 signalF.SetParameter(4, fgPedestalValue);
6751f762 861
862 Double_t signal=0.0, noise=0.0;
ee299369 863 for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
916f1e76 864 signal = signalF.Eval(iTime) ;
865
7643e728 866 // Next lines commeted for the moment but in principle it is not necessary to add
ff10f540 867 // extra noise since noise already added at the digits level.
7643e728 868
fe93d365 869 //According to Terry Awes, 13-Apr-2008
870 //add gaussian noise in quadrature to each sample
09974781 871 //Double_t noise = gRandom->Gaus(0.,fgFEENoise);
fe93d365 872 //signal = sqrt(signal*signal + noise*noise);
873
e2c2134b 874 // March 17,09 for fast fit simulations by Alexei Pavlinov.
4fe71e02 875 // Get from PHOS analysis. In some sense it is open questions.
6751f762 876 if(keyErr>0) {
877 noise = gRandom->Gaus(0.,fgFEENoise);
878 signal += noise;
879 }
880
ee299369 881 adcH[iTime] = static_cast<Int_t>(signal + 0.5) ;
b4133f05 882 if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits
883 adcH[iTime] = fgkRawSignalOverflow ;
ee299369 884 lowGain = kTRUE ;
885 }
886
887 signal /= fHighLowGainFactor;
888
889 adcL[iTime] = static_cast<Int_t>(signal + 0.5) ;
b4133f05 890 if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits
891 adcL[iTime] = fgkRawSignalOverflow ;
ee299369 892 }
893 return lowGain ;
894}
4fe71e02 895
829ba234 896//__________________________________________________________________
3a187ba0 897void AliEMCALRawUtils::CalculateChi2(const Double_t* t, const Double_t* y, const Int_t nPoints,
898const Double_t sig, const Double_t tau, const Double_t amp, const Double_t t0, Double_t &chi2)
899{
900 // Input:
901 // t[] - array of time bins
902 // y[] - array of amplitudes after pedestal subtractions;
903 // nPoints - number of points
904 // sig - error of amplitude measurement (one value for all channels)
905 // if sig<0 that mean sig=1.
906 // tau - filter time response (in timebin units)
907 // amp - amplitude at t0;
908 // t0 - time of max amplitude;
909 // Output:
910 // chi2 - chi2
911 // ndf = nPoints - 2 when tau fixed
912 // ndf = nPoints - 3 when tau free
913 static Double_t par[5]={0.0, 0.0, 0.0, 2.0, 0.0};
914
915 par[0] = amp;
916 par[1] = t0;
917 par[2] = tau;
918 // par[3]=n=2.; par[4]=ped=0.0
919
920 Double_t dy = 0.0, x = 0.0, f=0.0;
921 for(Int_t i=0; i<nPoints; i++){
922 x = t[i];
923 f = RawResponseFunction(&x, par);
924 dy = y[i] - f;
925 chi2 += dy*dy;
926 printf(" AliEMCALRawUtils::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, y[i], f, dy);
927 }
928 if(sig>0.0) chi2 /= (sig*sig);
929}
930
4fe71e02 931//__________________________________________________________________
932void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo)
933{
934 //Set fitting algorithm and initialize it if this same algorithm was not set before.
916f1e76 935 //printf("**** Set Algorithm , number %d ****\n",fitAlgo);
936
4fe71e02 937 if(fitAlgo == fFittingAlgorithm && fRawAnalyzer) {
938 //Do nothing, this same algorithm already set before.
939 //printf("**** Algorithm already set before, number %d, %s ****\n",fitAlgo, fRawAnalyzer->GetName());
940 return;
941 }
942 //Initialize the requested algorithm
943 if(fitAlgo != fFittingAlgorithm || !fRawAnalyzer) {
944 //printf("**** Init Algorithm , number %d ****\n",fitAlgo);
945
946 fFittingAlgorithm = fitAlgo;
947 if (fRawAnalyzer) delete fRawAnalyzer; // delete prev. analyzer if existed.
948
949 if (fitAlgo == kFastFit) {
950 fRawAnalyzer = new AliCaloRawAnalyzerFastFit();
951 }
952 else if (fitAlgo == kNeuralNet) {
953 fRawAnalyzer = new AliCaloRawAnalyzerNN();
954 }
955 else if (fitAlgo == kLMS) {
956 fRawAnalyzer = new AliCaloRawAnalyzerLMS();
957 }
958 else if (fitAlgo == kPeakFinder) {
959 fRawAnalyzer = new AliCaloRawAnalyzerPeakFinder();
960 }
961 else if (fitAlgo == kCrude) {
962 fRawAnalyzer = new AliCaloRawAnalyzerCrude();
963 }
964 else {
965 fRawAnalyzer = new AliCaloRawAnalyzer();
966 }
967 }
968
969}
970
971