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