Updated version.
[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 {
e853f058 111 AliDebug(1, 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
2cd0ffda 343 fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done, for the new methods
344 fRawAnalyzer->SetOverflowCut(fgkOverflowCut);
16605c06 345 fRawAnalyzer->SetAmpCut(fNoiseThreshold);
346 fRawAnalyzer->SetFitArrayCut(fNoiseThreshold);
347 fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later
ee299369 348
16605c06 349 // channel info parameters
829ba234 350 Int_t lowGain = 0;
e5bbbc4e 351 Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref.
ee299369 352
32cd4c24 353 // start loop over input stream
354 while (in.NextDDL()) {
916f1e76 355
356// if ( in.GetDDLNumber() != 0 && in.GetDDLNumber() != 2 ) continue;
357
32cd4c24 358 while (in.NextChannel()) {
7643e728 359
360 //Check if the signal is high or low gain and then do the fit,
16605c06 361 //if it is from TRU or LEDMon do not fit
7643e728 362 caloFlag = in.GetCaloFlag();
916f1e76 363// if (caloFlag != 0 && caloFlag != 1) continue;
364 if (caloFlag > 2) continue; // Work with ALTRO and FALTRO
365
366 //Do not fit bad channels of ALTRO
367 if(caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) {
5e3106bc 368 //printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow());
369 continue;
370 }
371
16605c06 372 vector<AliCaloBunchInfo> bunchlist;
32cd4c24 373 while (in.NextBunch()) {
16605c06 374 bunchlist.push_back( AliCaloBunchInfo(in.GetStartTimeBin(), in.GetBunchLength(), in.GetSignals() ) );
375 } // loop over bunches
7643e728 376
916f1e76 377
507751ce 378 if ( caloFlag < 2 ){ // ALTRO
916f1e76 379
507751ce 380 Float_t time = 0;
829ba234 381 Float_t amp = 0;
382 short timeEstimate = 0;
507751ce 383 Float_t ampEstimate = 0;
384 Bool_t fitDone = kFALSE;
2cd0ffda 385 Float_t chi2 = 0;
386 Int_t ndf = 0;
387
7683df1d 388 if ( fFittingAlgorithm == kFastFit || fFittingAlgorithm == kNeuralNet || fFittingAlgorithm == kLMS || fFittingAlgorithm == kPeakFinder || fFittingAlgorithm == kCrude) {
16605c06 389 // all functionality to determine amp and time etc is encapsulated inside the Evaluate call for these methods
390 AliCaloFitResults fitResults = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2());
391
829ba234 392 amp = fitResults.GetAmp();
393 time = fitResults.GetTime();
507751ce 394 timeEstimate = fitResults.GetMaxTimebin();
829ba234 395 ampEstimate = fitResults.GetMaxSig();
2cd0ffda 396 chi2 = fitResults.GetChi2();
397 ndf = fitResults.GetNdf();
507751ce 398 if (fitResults.GetStatus() == AliCaloFitResults::kFitPar) {
399 fitDone = kTRUE;
2cd0ffda 400 }
16605c06 401 }
402 else { // for the other methods we for now use the functionality of
403 // AliCaloRawAnalyzer as well, to select samples and prepare for fits,
404 // if it looks like there is something to fit
405
406 // parameters init.
507751ce 407 Float_t pedEstimate = 0;
16605c06 408 short maxADC = 0;
16605c06 409 Int_t first = 0;
410 Int_t last = 0;
411 Int_t bunchIndex = 0;
412 //
413 // The PreFitEvaluateSamples + later call to FitRaw will hopefully
414 // be replaced by a single Evaluate call or so soon, like for the other
415 // methods, but this should be good enough for evaluation of
416 // the methods for now (Jan. 2010)
417 //
418 int nsamples = fRawAnalyzer->PreFitEvaluateSamples( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2(), bunchIndex, ampEstimate, maxADC, timeEstimate, pedEstimate, first, last);
7643e728 419
1afbf4b1 420 if (ampEstimate >= fNoiseThreshold) { // something worth looking at
7643e728 421
f57baa2d 422 time = timeEstimate; // maxrev in AliCaloRawAnalyzer speak; comes with an offset w.r.t. real timebin
423 Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
16605c06 424 amp = ampEstimate;
425
2cd0ffda 426 if ( nsamples > 1 && maxADC<fgkOverflowCut ) { // possibly something to fit
427 FitRaw(first, last, amp, time, chi2, fitDone);
f57baa2d 428 time += timebinOffset;
507751ce 429 timeEstimate += timebinOffset;
2cd0ffda 430 ndf = nsamples - 2;
9f467289 431 }
16605c06 432
16605c06 433 } // ampEstimate check
434 } // method selection
507751ce 435
436 if ( fitDone ) { // brief sanity check of fit results
437 Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
438 Float_t timeDiff = time - timeEstimate;
439 if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) ) {
440 // AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations amp %f time %d", amp, time, ampEstimate, timeEstimate));
441
442 // for now just overwrite the fit results with the simple/initial estimate
829ba234 443 amp = ampEstimate;
444 time = timeEstimate;
507751ce 445 fitDone = kFALSE;
446 }
447 } // fitDone
16605c06 448
2cd0ffda 449 if (amp >= fNoiseThreshold) { // something to be stored
507751ce 450 if ( ! fitDone) { // smear ADC with +- 0.5 uniform (avoid discrete effects)
451 amp += (0.5 - gRandom->Rndm()); // Rndm generates a number in ]0,1]
452 }
453
829ba234 454 Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
455 lowGain = in.IsLowGain();
7643e728 456
16605c06 457 // go from time-bin units to physical time fgtimetrigger
458 time = time * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger?
e7acbc37 459 // subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0:
460 time -= in.GetL1Phase();
7643e728 461
462 AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
463 // printf("Added tower: SM %d, row %d, column %d, amp %3.2f\n",in.GetModule(), in.GetRow(), in.GetColumn(),amp);
2cd0ffda 464 AddDigit(digitsArr, id, lowGain, amp, time, chi2, ndf);
7643e728 465 }
466
916f1e76 467 }//ALTRO
46f1d25f 468 else if(fUseFALTRO)
916f1e76 469 {// Fake ALTRO
470 // if (maxTimeBin && gSig->GetN() > maxTimeBin + 10) gSig->Set(maxTimeBin + 10); // set actual max size of TGraph
916f1e76 471 Int_t hwAdd = in.GetHWAddress();
472 UShort_t iRCU = in.GetDDLNumber() % 2; // 0/1
473 UShort_t iBranch = ( hwAdd >> 11 ) & 0x1; // 0/1
474
475 // Now find TRU number
476 Int_t itru = 3 * in.GetModule() + ( (iRCU << 1) | iBranch ) - 1;
477
478 AliDebug(1,Form("Found TRG digit in TRU: %2d ADC: %2d",itru,in.GetColumn()));
479
480 Int_t idtrg;
481
482 Bool_t isOK = fGeom->GetAbsFastORIndexFromTRU(itru, in.GetColumn(), idtrg);
483
484 Int_t timeSamples[256]; for (Int_t j=0;j<256;j++) timeSamples[j] = 0;
485 Int_t nSamples = 0;
486
487 for (std::vector<AliCaloBunchInfo>::iterator itVectorData = bunchlist.begin(); itVectorData != bunchlist.end(); itVectorData++)
488 {
489 AliCaloBunchInfo bunch = *(itVectorData);
490
491 const UShort_t* sig = bunch.GetData();
492 Int_t startBin = bunch.GetStartBin();
493
494 for (Int_t iS = 0; iS < bunch.GetLength(); iS++)
495 {
496 Int_t time = startBin--;
497 Int_t amp = sig[iS];
498
499 if ( amp ) timeSamples[nSamples++] = ( ( time << 12 ) & 0xFF000 ) | ( amp & 0xFFF );
500 }
501 }
502
503 if (nSamples && isOK) AddDigit(digitsTRG, idtrg, timeSamples, nSamples);
504 }//Fake ALTRO
32cd4c24 505 } // end while over channel
506 } //end while over DDL's, of input stream
16605c06 507
f4a4dc82 508 TrimDigits(digitsArr);
509
ee299369 510 return ;
511}
512
513//____________________________________________________________________________
916f1e76 514void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t timeSamples[], Int_t nSamples)
515{
f4a4dc82 516 //Add raw sample to raw digit
517 new((*digitsArr)[digitsArr->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
518
519 // Int_t idx = digitsArr->GetEntriesFast()-1;
520 // AliEMCALRawDigit* d = (AliEMCALRawDigit*)digitsArr->At(idx);
916f1e76 521}
522
523//____________________________________________________________________________
2cd0ffda 524void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf) {
82cbdfca 525 //
526 // Add a new digit.
527 // This routine checks whether a digit exists already for this tower
528 // and then decides whether to use the high or low gain info
529 //
530 // Called by Raw2Digits
531
532 AliEMCALDigit *digit = 0, *tmpdigit = 0;
82cbdfca 533 TIter nextdigit(digitsArr);
534 while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) {
f4a4dc82 535 if (tmpdigit->GetId() == id) digit = tmpdigit;
82cbdfca 536 }
537
538 if (!digit) { // no digit existed for this tower; create one
f4a4dc82 539 Int_t type = AliEMCALDigit::kHG; // use enum in AliEMCALDigit
540 if (lowGain) {
541 amp *= fHighLowGainFactor;
542 type = AliEMCALDigit::kLGnoHG;
543 }
544 Int_t idigit = digitsArr->GetEntries();
2cd0ffda 545 new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit, chi2, ndf);
f4a4dc82 546 AliDebug(2,Form("Add digit Id %d for the first time, type %d", id, type));
547 }//digit added first time
82cbdfca 548 else { // a digit already exists, check range
f4a4dc82 549 // (use high gain if signal < cut value, otherwise low gain)
550 if (lowGain) { // new digit is low gain
551 if (digit->GetAmplitude() > fgkOverflowCut) { // use if previously stored (HG) digit is out of range
552 digit->SetAmplitude(fHighLowGainFactor * amp);
553 digit->SetTime(time);
554 digit->SetType(AliEMCALDigit::kLG);
555 AliDebug(2,Form("Add LG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType()));
556 }
557 }//new low gain digit
558 else { // new digit is high gain
559 if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range
560 digit->SetAmplitude(amp);
561 digit->SetTime(time);
562 digit->SetType(AliEMCALDigit::kHG);
563 AliDebug(2,Form("Add HG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType()));
564 }
565 else { // HG out of range, just change flag value to show that HG did exist
566 digit->SetType(AliEMCALDigit::kLG);
567 AliDebug(2,Form("Change LG digit to HG, ID %d, type %d", digit->GetId(), digit->GetType()));
568 }
569 }//new high gain digit
570 }//digit existed replace it
571
572}
573
574//____________________________________________________________________________
575void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr)
576{
577 // Remove digits with only low gain and large time
578
579 AliEMCALDigit *digit = 0;
580 Int_t n = 0;
581 Int_t nDigits = digitsArr->GetEntriesFast();
582 TIter nextdigit(digitsArr);
583 while ((digit = (AliEMCALDigit*) nextdigit())) {
584
585 //Check if only LG existed, remove if so
586 if (digit->GetType() == AliEMCALDigit::kLGnoHG) {
587 AliDebug(1,Form("Remove digit with id %d, LGnoHG",digit->GetId()));
588 digitsArr->Remove(digit);
82cbdfca 589 }
2cd0ffda 590 //Check if time is too large or too small, remove if so
f4a4dc82 591 else if(fTimeMin > digit->GetTime() || fTimeMax < digit->GetTime()) {
592 digitsArr->Remove(digit);
593 AliDebug(1,Form("Remove digit with id %d, Bad Time %e",digit->GetId(), digit->GetTime()));
82cbdfca 594 }
2cd0ffda 595 // Check if Chi2 is undefined
596 else if (0 > digit->GetChi2()) {
597 digitsArr->Remove(digit);
598 AliDebug(1,Form("Remove digit with id %d, Bad Chi2 %e",digit->GetId(), digit->GetChi2()));
599 }
f4a4dc82 600 //Good digit, just reassign the index of the digit in case there was a previous removal
601 else {
602 digit->SetIndexInList(n);
603 n++;
604 }
605 }//while
606
607 digitsArr->Compress();
608 AliDebug(1,Form("N Digits before trimming : %d; after array compression %d",nDigits,digitsArr->GetEntriesFast()));
609
82cbdfca 610}
f4a4dc82 611
82cbdfca 612//____________________________________________________________________________
2cd0ffda 613void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
16605c06 614{ // Fits the raw signal time distribution
615
616 //--------------------------------------------------
617 //Do the fit, different fitting algorithms available
618 //--------------------------------------------------
619 int nsamples = lastTimeBin - firstTimeBin + 1;
507751ce 620 fitDone = kFALSE;
ee299369 621
16605c06 622 switch(fFittingAlgorithm) {
623 case kStandard:
624 {
7683df1d 625 if (nsamples < 3) { return; } // nothing much to fit
16605c06 626 //printf("Standard fitter \n");
7683df1d 627
16605c06 628 // Create Graph to hold data we will fit
7683df1d 629 TGraph *gSig = new TGraph( nsamples);
630 for (int i=0; i<nsamples; i++) {
631 Int_t timebin = firstTimeBin + i;
f57baa2d 632 gSig->SetPoint(i, timebin, fRawAnalyzer->GetReversed(timebin));
7683df1d 633 }
634
16605c06 635 TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
636 signalF->SetParameters(10.,5.,fTau,fOrder,0.); //set all defaults once, just to be safe
637 signalF->SetParNames("amp","t0","tau","N","ped");
638 signalF->FixParameter(2,fTau); // tau in units of time bin
639 signalF->FixParameter(3,fOrder); // order
640 signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here
641 signalF->SetParameter(1, time);
642 signalF->SetParameter(0, amp);
507751ce 643 // set rather loose parameter limits
644 signalF->SetParLimits(0, 0.5*amp, 2*amp );
645 signalF->SetParLimits(1, time - 4, time + 4);
646
647 try {
648 gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
649 // assign fit results
829ba234 650 amp = signalF->GetParameter(0);
507751ce 651 time = signalF->GetParameter(1);
2cd0ffda 652 chi2 = signalF->GetChisquare();
507751ce 653 fitDone = kTRUE;
654 }
655 catch (const std::exception & e) {
656 AliError( Form("TGraph Fit exception %s", e.what()) );
657 // stay with default amp and time in case of exception, i.e. no special action required
658 fitDone = kFALSE;
659 }
16605c06 660 delete signalF;
661
16605c06 662 //printf("Std : Amp %f, time %g\n",amp, time);
7683df1d 663 delete gSig; // delete TGraph
16605c06 664
665 break;
666 }//kStandard Fitter
667 //----------------------------
7683df1d 668 case kLogFit:
16605c06 669 {
7683df1d 670 if (nsamples < 3) { return; } // nothing much to fit
671 //printf("LogFit \n");
672
673 // Create Graph to hold data we will fit
674 TGraph *gSigLog = new TGraph( nsamples);
675 for (int i=0; i<nsamples; i++) {
676 Int_t timebin = firstTimeBin + i;
677 gSigLog->SetPoint(timebin, timebin, TMath::Log(fRawAnalyzer->GetReversed(timebin) ) );
7643e728 678 }
7683df1d 679
680 TF1 * signalFLog = new TF1("signalLog", RawResponseFunctionLog, 0, GetRawFormatTimeBins(), 5);
681 signalFLog->SetParameters(2.3, 5.,fTau,fOrder,0.); //set all defaults once, just to be safe
682 signalFLog->SetParNames("amplog","t0","tau","N","ped");
683 signalFLog->FixParameter(2,fTau); // tau in units of time bin
684 signalFLog->FixParameter(3,fOrder); // order
685 signalFLog->FixParameter(4, 0); // pedestal should be subtracted when we get here
686 signalFLog->SetParameter(1, time);
687 if (amp>=1) {
688 signalFLog->SetParameter(0, TMath::Log(amp));
16605c06 689 }
7683df1d 690
691 gSigLog->Fit(signalFLog, "QROW"); // Note option 'W': equal errors on all points
692
693 // assign fit results
694 Double_t amplog = signalFLog->GetParameter(0); //Not Amp, but Log of Amp
695 amp = TMath::Exp(amplog);
696 time = signalFLog->GetParameter(1);
507751ce 697 fitDone = kTRUE;
7683df1d 698
699 delete signalFLog;
700 //printf("LogFit: Amp %f, time %g\n",amp, time);
701 delete gSigLog;
16605c06 702 break;
7683df1d 703 } //kLogFit
704 //----------------------------
705
16605c06 706 //----------------------------
707 }//switch fitting algorithms
fb070798 708
16605c06 709 return;
710}
8cb998bd 711
16605c06 712//__________________________________________________________________
713void AliEMCALRawUtils::FitParabola(const TGraph *gSig, Float_t & amp) const
714{
715 //BEG YS alternative methods to calculate the amplitude
716 Double_t * ymx = gSig->GetX() ;
717 Double_t * ymy = gSig->GetY() ;
718 const Int_t kN = 3 ;
719 Double_t ymMaxX[kN] = {0., 0., 0.} ;
720 Double_t ymMaxY[kN] = {0., 0., 0.} ;
721 Double_t ymax = 0. ;
722 // find the maximum amplitude
723 Int_t ymiMax = 0 ;
724 for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
725 if (ymy[ymi] > ymMaxY[0] ) {
726 ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
727 ymMaxX[0] = ymx[ymi] ;
728 ymiMax = ymi ;
729 }
730 }
731 // find the maximum by fitting a parabola through the max and the two adjacent samples
732 if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
733 ymMaxY[1] = ymy[ymiMax+1] ;
734 ymMaxY[2] = ymy[ymiMax-1] ;
735 ymMaxX[1] = ymx[ymiMax+1] ;
736 ymMaxX[2] = ymx[ymiMax-1] ;
737 if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
738 //fit a parabola through the 3 points y= a+bx+x*x*x
739 Double_t sy = 0 ;
740 Double_t sx = 0 ;
741 Double_t sx2 = 0 ;
742 Double_t sx3 = 0 ;
743 Double_t sx4 = 0 ;
744 Double_t sxy = 0 ;
745 Double_t sx2y = 0 ;
746 for (Int_t i = 0; i < kN ; i++) {
747 sy += ymMaxY[i] ;
748 sx += ymMaxX[i] ;
749 sx2 += ymMaxX[i]*ymMaxX[i] ;
750 sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
751 sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
752 sxy += ymMaxX[i]*ymMaxY[i] ;
753 sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ;
754 }
755 Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2);
756 Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
757 Double_t c = cN / cD ;
758 Double_t b = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
759 Double_t a = (sy-b*sx-c*sx2)/kN ;
760 Double_t xmax = -b/(2*c) ;
761 ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
029fe7a2 762 amp = ymax;
16605c06 763 }
764 }
765
766 Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ;
767 if (diff > 0.1)
768 amp = ymMaxY[0] ;
769 //printf("Yves : Amp %f, time %g\n",amp, time);
770 //END YS
ee299369 771 return;
772}
16605c06 773
ee299369 774//__________________________________________________________________
775Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
776{
8cb998bd 777 // Matches version used in 2007 beam test
778 //
ee299369 779 // Shape of the electronics raw reponse:
780 // It is a semi-gaussian, 2nd order Gamma function of the general form
781 //
7643e728 782 // xx = (t - t0 + tau) / tau [xx is just a convenient help variable]
783 // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0
784 // F = 0 for xx < 0
ee299369 785 //
786 // parameters:
8cb998bd 787 // A: par[0] // Amplitude = peak value
788 // t0: par[1]
789 // tau: par[2]
790 // N: par[3]
791 // ped: par[4]
ee299369 792 //
793 Double_t signal ;
8cb998bd 794 Double_t tau =par[2];
e5bbbc4e 795 Double_t n =par[3];
8cb998bd 796 Double_t ped = par[4];
797 Double_t xx = ( x[0] - par[1] + tau ) / tau ;
ee299369 798
5a056daa 799 if (xx <= 0)
8cb998bd 800 signal = ped ;
ee299369 801 else {
e5bbbc4e 802 signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
ee299369 803 }
804 return signal ;
805}
806
807//__________________________________________________________________
7683df1d 808Double_t AliEMCALRawUtils::RawResponseFunctionLog(Double_t *x, Double_t *par)
809{
810 // Matches version used in 2007 beam test
811 //
812 // Shape of the electronics raw reponse:
813 // It is a semi-gaussian, 2nd order Gamma function of the general form
814 //
815 // xx = (t - t0 + tau) / tau [xx is just a convenient help variable]
816 // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0
817 // F = 0 for xx < 0
818 //
819 // parameters:
820 // Log[A]: par[0] // Amplitude = peak value
821 // t0: par[1]
822 // tau: par[2]
823 // N: par[3]
824 // ped: par[4]
825 //
826 Double_t signal ;
827 Double_t tau =par[2];
828 Double_t n =par[3];
829 //Double_t ped = par[4]; // not used
830 Double_t xx = ( x[0] - par[1] + tau ) / tau ;
831
832 if (xx < 0)
833 signal = par[0] - n*TMath::Log(TMath::Abs(xx)) + n * (1 - xx ) ;
834 else {
835 signal = par[0] + n*TMath::Log(xx) + n * (1 - xx ) ;
836 }
837 return signal ;
838}
839
840//__________________________________________________________________
6751f762 841Bool_t AliEMCALRawUtils::RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL, const Int_t keyErr) const
ee299369 842{
843 // for a start time dtime and an amplitude damp given by digit,
844 // calculates the raw sampled response AliEMCAL::RawResponseFunction
845
ee299369 846 Bool_t lowGain = kFALSE ;
847
48a56166 848 // A: par[0] // Amplitude = peak value
849 // t0: par[1]
850 // tau: par[2]
851 // N: par[3]
852 // ped: par[4]
853
56e13066 854 TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
48a56166 855 signalF.SetParameter(0, damp) ;
56e13066 856 signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ;
b4133f05 857 signalF.SetParameter(2, fTau) ;
858 signalF.SetParameter(3, fOrder);
fe93d365 859 signalF.SetParameter(4, fgPedestalValue);
6751f762 860
861 Double_t signal=0.0, noise=0.0;
ee299369 862 for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
916f1e76 863 signal = signalF.Eval(iTime) ;
864
7643e728 865 // Next lines commeted for the moment but in principle it is not necessary to add
ff10f540 866 // extra noise since noise already added at the digits level.
7643e728 867
fe93d365 868 //According to Terry Awes, 13-Apr-2008
869 //add gaussian noise in quadrature to each sample
09974781 870 //Double_t noise = gRandom->Gaus(0.,fgFEENoise);
fe93d365 871 //signal = sqrt(signal*signal + noise*noise);
872
e2c2134b 873 // March 17,09 for fast fit simulations by Alexei Pavlinov.
4fe71e02 874 // Get from PHOS analysis. In some sense it is open questions.
6751f762 875 if(keyErr>0) {
876 noise = gRandom->Gaus(0.,fgFEENoise);
877 signal += noise;
878 }
879
ee299369 880 adcH[iTime] = static_cast<Int_t>(signal + 0.5) ;
b4133f05 881 if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits
882 adcH[iTime] = fgkRawSignalOverflow ;
ee299369 883 lowGain = kTRUE ;
884 }
885
886 signal /= fHighLowGainFactor;
887
888 adcL[iTime] = static_cast<Int_t>(signal + 0.5) ;
b4133f05 889 if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits
890 adcL[iTime] = fgkRawSignalOverflow ;
ee299369 891 }
892 return lowGain ;
893}
4fe71e02 894
829ba234 895//__________________________________________________________________
3a187ba0 896void AliEMCALRawUtils::CalculateChi2(const Double_t* t, const Double_t* y, const Int_t nPoints,
897const Double_t sig, const Double_t tau, const Double_t amp, const Double_t t0, Double_t &chi2)
898{
899 // Input:
900 // t[] - array of time bins
901 // y[] - array of amplitudes after pedestal subtractions;
902 // nPoints - number of points
903 // sig - error of amplitude measurement (one value for all channels)
904 // if sig<0 that mean sig=1.
905 // tau - filter time response (in timebin units)
906 // amp - amplitude at t0;
907 // t0 - time of max amplitude;
908 // Output:
909 // chi2 - chi2
910 // ndf = nPoints - 2 when tau fixed
911 // ndf = nPoints - 3 when tau free
912 static Double_t par[5]={0.0, 0.0, 0.0, 2.0, 0.0};
913
914 par[0] = amp;
915 par[1] = t0;
916 par[2] = tau;
917 // par[3]=n=2.; par[4]=ped=0.0
918
919 Double_t dy = 0.0, x = 0.0, f=0.0;
920 for(Int_t i=0; i<nPoints; i++){
921 x = t[i];
922 f = RawResponseFunction(&x, par);
923 dy = y[i] - f;
924 chi2 += dy*dy;
925 printf(" AliEMCALRawUtils::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, y[i], f, dy);
926 }
927 if(sig>0.0) chi2 /= (sig*sig);
928}
929
4fe71e02 930//__________________________________________________________________
931void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo)
932{
933 //Set fitting algorithm and initialize it if this same algorithm was not set before.
916f1e76 934 //printf("**** Set Algorithm , number %d ****\n",fitAlgo);
935
4fe71e02 936 if(fitAlgo == fFittingAlgorithm && fRawAnalyzer) {
937 //Do nothing, this same algorithm already set before.
938 //printf("**** Algorithm already set before, number %d, %s ****\n",fitAlgo, fRawAnalyzer->GetName());
939 return;
940 }
941 //Initialize the requested algorithm
942 if(fitAlgo != fFittingAlgorithm || !fRawAnalyzer) {
943 //printf("**** Init Algorithm , number %d ****\n",fitAlgo);
944
945 fFittingAlgorithm = fitAlgo;
946 if (fRawAnalyzer) delete fRawAnalyzer; // delete prev. analyzer if existed.
947
948 if (fitAlgo == kFastFit) {
949 fRawAnalyzer = new AliCaloRawAnalyzerFastFit();
950 }
951 else if (fitAlgo == kNeuralNet) {
952 fRawAnalyzer = new AliCaloRawAnalyzerNN();
953 }
954 else if (fitAlgo == kLMS) {
955 fRawAnalyzer = new AliCaloRawAnalyzerLMS();
956 }
957 else if (fitAlgo == kPeakFinder) {
958 fRawAnalyzer = new AliCaloRawAnalyzerPeakFinder();
959 }
960 else if (fitAlgo == kCrude) {
961 fRawAnalyzer = new AliCaloRawAnalyzerCrude();
962 }
963 else {
964 fRawAnalyzer = new AliCaloRawAnalyzer();
965 }
966 }
967
968}
969
970