Line data Source code
1 : // -*- mode: c++ -*-
2 : /**************************************************************************
3 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 : * *
5 : * Author: The ALICE Off-line Project. *
6 : * Contributors are mentioned in the code where appropriate. *
7 : * *
8 : * Permission to use, copy, modify and distribute this software and its *
9 : * documentation strictly for non-commercial purposes is hereby granted *
10 : * without fee, provided that the above copyright notice appears in all *
11 : * copies and that both the copyright notice and this permission notice *
12 : * appear in the supporting documentation. The authors make no claims *
13 : * about the suitability of this software for any purpose. It is *
14 : * provided "as is" without express or implied warranty. *
15 : **************************************************************************/
16 :
17 : //#include "AliCDBManager.h"
18 : #include "AliEMCALRawUtils.h"
19 : #include "AliRun.h"
20 : #include "AliRunLoader.h"
21 : #include "AliAltroBuffer.h"
22 : #include "AliRawReader.h"
23 : #include "AliCaloRawStreamV3.h"
24 : #include "AliDAQ.h"
25 : #include "AliEMCALRecParam.h"
26 : #include "AliEMCALLoader.h"
27 : #include "AliEMCALGeometry.h"
28 : #include "AliEMCALDigit.h"
29 : #include "AliEMCALRawDigit.h"
30 : #include "AliEMCAL.h"
31 : #include "AliCaloCalibPedestal.h"
32 : #include "AliCaloBunchInfo.h"
33 : #include "AliCaloFitResults.h"
34 : #include "AliEMCALTriggerRawDigitMaker.h"
35 : #include "AliEMCALTriggerSTURawStream.h"
36 : #include "AliEMCALTriggerData.h"
37 : #include "AliCaloConstants.h"
38 : #include "AliCaloRawAnalyzer.h"
39 : #include "AliCaloRawAnalyzerFactory.h"
40 : #include "AliEMCALRawResponse.h"
41 :
42 : using namespace CALO;
43 : using namespace EMCAL;
44 :
45 : using std::vector;
46 :
47 : /// \cond CLASSIMP
48 42 : ClassImp(AliEMCALRawUtils) ;
49 : /// \endcond
50 :
51 : ///
52 : /// Constructor. Set up fitting algorightm, geometry
53 : /// and default parameter values.
54 : ///
55 12 : AliEMCALRawUtils::AliEMCALRawUtils( Algo::fitAlgorithm fitAlgo) : fNoiseThreshold(3),
56 6 : fNPedSamples(4),
57 6 : fGeom(0),
58 6 : fOption(""),
59 6 : fRemoveBadChannels(kFALSE),
60 6 : fFittingAlgorithm(fitAlgo),
61 6 : fTimeMin(-1.),
62 6 : fTimeMax(1.),
63 6 : fUseFALTRO(kTRUE),
64 6 : fUseL1Phase(kTRUE),
65 6 : fRawAnalyzer(0),
66 6 : fTriggerRawDigitMaker(0x0)
67 30 : {
68 6 : SetFittingAlgorithm(fitAlgo);
69 :
70 6 : const TObjArray* maps = AliEMCALRecParam::GetMappings();
71 :
72 6 : if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
73 :
74 60 : for(Int_t i = 0; i < 4; i++)
75 : {
76 48 : fMapping[i] = (AliAltroMapping*)maps->At(i);
77 : }
78 :
79 6 : AliRunLoader *rl = AliRunLoader::Instance();
80 18 : if (rl && rl->GetAliRun())
81 : {
82 25 : AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
83 5 : if(emcal)
84 : {
85 10 : fGeom = emcal->GetGeometry();
86 5 : }
87 : else
88 : {
89 0 : AliDebug(1, Form("Using default geometry in raw reco"));
90 0 : fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
91 : }
92 5 : }
93 : else
94 : {
95 5 : AliDebug(1, Form("Using default geometry in raw reco"));
96 2 : fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
97 : }
98 :
99 6 : if(!fGeom) AliFatal(Form("Could not get geometry!"));
100 18 : fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();
101 12 : }
102 :
103 : ///
104 : /// Destructor.
105 : ///
106 : AliEMCALRawUtils::~AliEMCALRawUtils()
107 28 : {
108 12 : delete fRawAnalyzer;
109 12 : delete fTriggerRawDigitMaker;
110 14 : }
111 :
112 : ///
113 : /// Convert digits of the current event to raw data.
114 : ///
115 : void AliEMCALRawUtils::Digits2Raw()
116 : {
117 8 : AliRunLoader *rl = AliRunLoader::Instance();
118 :
119 12 : AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
120 4 : loader->LoadDigits("EMCAL");
121 4 : loader->GetEvent();
122 :
123 4 : TClonesArray* digits = loader->Digits() ;
124 :
125 4 : if (!digits)
126 : {
127 0 : AliWarning("No digits found !");
128 0 : return;
129 : }
130 :
131 : static const Int_t nDDL = 20*2; // 20 SM for EMCal + DCal hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
132 :
133 4 : AliAltroBuffer* buffers[nDDL];
134 328 : for (Int_t i=0; i < nDDL; i++)
135 160 : buffers[i] = 0;
136 :
137 4 : TArrayI adcValuesLow ( TIMEBINS );
138 4 : TArrayI adcValuesHigh( TIMEBINS );
139 :
140 : // Loop over digits (assume ordered digits)
141 186 : for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++)
142 : {
143 232 : AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
144 :
145 58 : if(!digit)
146 : {
147 0 : AliFatal("NULL Digit");
148 : }
149 : else
150 : {
151 58 : if (digit->GetAmplitude() < AliEMCALRawResponse::GetRawFormatThreshold() )
152 : {
153 0 : continue;
154 : }
155 :
156 : // Get cell indices
157 58 : Int_t nSM = 0;
158 58 : Int_t nIphi = 0;
159 58 : Int_t nIeta = 0;
160 58 : Int_t iphi = 0;
161 58 : Int_t ieta = 0;
162 58 : Int_t nModule = 0;
163 :
164 58 : fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
165 58 : fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
166 :
167 : //----------------------------------------------------------------------
168 : //
169 : // Online mapping and numbering is the same for EMCal and DCal SMs but:
170 : // - DCal odd SM (13,15,17) has online cols: 16-47; offline cols 0-31.
171 : // - Even DCal SMs have the same numbering online and offline 0-31.
172 : // - DCal 1/3 SM (18,19), online rows 16-23; offline rows 0-7
173 : //
174 : // In the next lines shift the online cols or rows depending on the
175 : // SM to match the offline mapping.
176 : // Apply the shifts (inverse to those in Raw2Digits):
177 :
178 58 : fGeom->ShiftOfflineToOnlineCellIndexes(nSM, iphi, ieta);
179 :
180 : //
181 : //----------------------------------------------------------------------
182 :
183 : // Check which is the RCU, 0 or 1, of the cell.
184 : // *** RCU corresponds to DDL in a SuperModule ***
185 : Int_t iRCU = -111;
186 64 : if ( 0<=iphi&&iphi< 8 ) iRCU=0; // first cable row
187 68 : else if ( 8<=iphi&&iphi<16 && 0<=ieta&&ieta<24 ) iRCU=0; // first half;
188 41 : else if ( 8<=iphi&&iphi<16 && 24<=ieta&&ieta<48 ) iRCU=1; // second half;
189 62 : else if ( 16<=iphi&&iphi<24 ) iRCU=1; // third cable row
190 :
191 106 : if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
192 :
193 58 : if (iRCU<0)
194 0 : Fatal("Digits2Raw()","Non-existent RCU/DDL number: %d", iRCU);
195 :
196 : // Which DDL?
197 58 : Int_t iDDL = NRCUSPERMODULE*nSM + iRCU;
198 :
199 58 : if (iDDL < 0 || iDDL >= nDDL)
200 : {
201 0 : Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
202 : }
203 : else
204 : {
205 58 : if (buffers[iDDL] == 0)
206 : {
207 : // open new file and write dummy header
208 40 : TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
209 : //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C
210 20 : Int_t iRCUside=iRCU+(nSM%2)*2;
211 : //iRCU=0 and even (0) SM -> RCU0A.data 0
212 : //iRCU=1 and even (0) SM -> RCU1A.data 1
213 : //iRCU=0 and odd (1) SM -> RCU0C.data 2
214 : //iRCU=1 and odd (1) SM -> RCU1C.data 3
215 80 : buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]);
216 20 : buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy;
217 20 : }
218 :
219 : // out of time range signal (?)
220 58 : if (digit->GetTimeR() > TIMEBINMAX )
221 : {
222 0 : AliInfo("Signal is out of time range.\n");
223 0 : buffers[iDDL]->FillBuffer((Int_t)digit->GetAmplitude());
224 0 : buffers[iDDL]->FillBuffer( TIMEBINS ); // time bin
225 0 : buffers[iDDL]->FillBuffer(3); // bunch length
226 0 : buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer
227 : // calculate the time response function
228 : }
229 : else
230 : {
231 174 : Bool_t lowgain = AliEMCALRawResponse::RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(),
232 58 : adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ;
233 :
234 116 : if (lowgain)
235 58 : buffers[iDDL]->WriteChannel(ieta, iphi, 0, TIMEBINS, adcValuesLow .GetArray(), AliEMCALRawResponse::GetRawFormatThreshold() );
236 : else
237 58 : buffers[iDDL]->WriteChannel(ieta, iphi, 1, TIMEBINS, adcValuesHigh.GetArray(), AliEMCALRawResponse::GetRawFormatThreshold() );
238 : }
239 : }// iDDL under the limits
240 58 : }// Digit exists
241 58 : }// Digit loop
242 :
243 : // Write headers and close files
244 328 : for (Int_t i=0; i < nDDL; i++)
245 : {
246 160 : if (buffers[i])
247 : {
248 20 : buffers[i]->Flush();
249 20 : buffers[i]->WriteDataHeader(kFALSE, kFALSE);
250 40 : delete buffers[i];
251 : }
252 : }
253 :
254 4 : loader->UnloadDigits();
255 8 : }
256 :
257 : ///
258 : /// Create the digit from raw fit and add it to the list of digits
259 : ///
260 : void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain,
261 : Float_t amp, Float_t time, Float_t chi2, Int_t ndf)
262 : {
263 : AliEMCALDigit *digit = 0, *tmpdigit = 0;
264 116 : TIter nextdigit(digitsArr);
265 :
266 1676 : while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit()))
267 : {
268 483 : if (tmpdigit->GetId() == id) digit = tmpdigit;
269 : }
270 :
271 : // No digit existed for this tower; create one.
272 58 : if (!digit)
273 : {
274 : Int_t type = AliEMCALDigit::kHG; // use enum in AliEMCALDigit
275 :
276 57 : if (lowGain)
277 : {
278 0 : amp *= HGLGFACTOR;
279 : type = AliEMCALDigit::kLGnoHG;
280 0 : }
281 :
282 57 : Int_t idigit = digitsArr->GetEntries();
283 :
284 171 : new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit, chi2, ndf);
285 :
286 285 : AliDebug(2,Form("Add digit Id %d for the first time, type %d", id, type));
287 57 : }// digit added first time.
288 :
289 : // A digit already exists, check range
290 : // (use high gain if signal < cut value, otherwise low gain)
291 : else
292 : {
293 1 : if (lowGain)
294 : {
295 : // New digit is low gain
296 0 : if (digit->GetAmplitude() > OVERFLOWCUT )
297 : {
298 : // Use if previously stored (HG) digit is out of range
299 0 : digit->SetAmplitude( HGLGFACTOR * amp);
300 0 : digit->SetTime(time);
301 0 : digit->SetType(AliEMCALDigit::kLG);
302 :
303 0 : AliDebug(2,Form("Add LG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType()));
304 : }
305 : } // New low gain digit
306 : else
307 : {
308 : // New digit is high gain
309 1 : if (amp < OVERFLOWCUT )
310 : {
311 : // New digit is high gain; use if not out of range
312 1 : digit->SetAmplitude(amp);
313 1 : digit->SetTime(time);
314 1 : digit->SetType(AliEMCALDigit::kHG);
315 :
316 5 : AliDebug(2,Form("Add HG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType()));
317 : }
318 : else
319 : {
320 : // HG out of range, just change flag value to show that HG did exist
321 0 : digit->SetType(AliEMCALDigit::kLG);
322 :
323 0 : AliDebug(2,Form("Change LG digit to HG, ID %d, type %d", digit->GetId(), digit->GetType()));
324 : }
325 : } // New high gain digit
326 : }// Digit existed replace it
327 58 : }
328 :
329 : ///
330 : /// Conversion of raw data to digits.
331 : ///
332 : void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap,
333 : TClonesArray *digitsTRG, TClonesArray *trgData)
334 : {
335 12 : if ( digitsArr) digitsArr->Clear("C");
336 :
337 4 : if (!digitsArr) { Error("Raw2Digits", "no digits found !") ; return ; }
338 :
339 4 : if (!reader) { Error("Raw2Digits", "no raw reader found !"); return ; }
340 :
341 4 : AliEMCALTriggerSTURawStream inSTU(reader);
342 :
343 12 : AliCaloRawStreamV3 in(reader,"EMCAL",fMapping);
344 :
345 4 : reader->Select("EMCAL",0,AliDAQ::GetFirstSTUDDL()-1);
346 :
347 4 : fTriggerRawDigitMaker->Reset();
348 4 : fTriggerRawDigitMaker->SetIO(reader, in, inSTU, digitsTRG, trgData);
349 :
350 4 : fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later
351 :
352 : Int_t lowGain = 0;
353 : Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref.
354 :
355 : Float_t bcTimePhaseCorr = 0; // for BC-based L1 phase correction, Run1 data
356 8 : Int_t bcMod4 = (reader->GetBCID() % 4); // LHC uses 40 MHz, EMCal uses 10 MHz clock
357 :
358 : //AliCDBManager* man = AliCDBManager::Instance();
359 : //Int_t runNumber = man->GetRun();
360 :
361 4 : Int_t runNumber = reader->GetRunNumber();
362 :
363 : // Apply this shift for Run1 data
364 4 : if ((runNumber > 130850 && runNumber < 200000) && (bcMod4==0 || bcMod4==1))
365 0 : bcTimePhaseCorr = -1e-7; // subtract 100 ns for certain BC values
366 :
367 48 : while (in.NextDDL())
368 : {
369 156 : while (in.NextChannel())
370 : {
371 58 : caloFlag = in.GetCaloFlag();
372 :
373 58 : if ( caloFlag > 2 ) continue; // Work with ALTRO and FALTRO
374 :
375 :
376 58 : Int_t sm = in.GetModule() ;
377 58 : Int_t row = in.GetRow () ;
378 58 : Int_t column = in.GetColumn() ;
379 :
380 : //----------------------------------------------------------------------
381 : //
382 : // Online mapping and numbering is the same for EMCal and DCal SMs but:
383 : // - DCal odd SM (13,15,17) has online cols: 16-47; offline cols 0-31.
384 : // - Even DCal SMs have the same numbering online and offline 0-31.
385 : // - DCal 1/3 SM (18,19), online rows 16-23; offline rows 0-7
386 : //
387 : // In the next lines shift the online cols or rows depending on the
388 : // SM to match the offline mapping.
389 : //
390 :
391 58 : fGeom->ShiftOnlineToOfflineCellIndexes(sm, row, column);
392 :
393 : //
394 : //---------------------------------------------------------------------
395 :
396 116 : if ( caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(sm, column, row) )
397 : {
398 0 : continue;
399 : }
400 :
401 58 : vector<AliCaloBunchInfo> bunchlist;
402 :
403 348 : while (in.NextBunch())
404 : {
405 174 : bunchlist.push_back( AliCaloBunchInfo(in.GetStartTimeBin(), in.GetBunchLength(), in.GetSignals() ) );
406 : }
407 :
408 58 : if (bunchlist.size() == 0) continue;
409 :
410 58 : if ( caloFlag < 2 )
411 : {
412 : // ALTRO
413 :
414 58 : Int_t id = fGeom->GetAbsCellIdFromCellIndexes(sm, row, column) ;
415 :
416 58 : lowGain = in.IsLowGain();
417 :
418 116 : if(fUseL1Phase)
419 174 : fRawAnalyzer->SetL1Phase( in.GetL1Phase() );
420 : else
421 0 : fRawAnalyzer->SetL1Phase( 0 );
422 :
423 58 : AliCaloFitResults res = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2());
424 :
425 58 : if(res.GetAmp() >= fNoiseThreshold )
426 : {
427 58 : AddDigit(digitsArr, id, lowGain, res.GetAmp(), res.GetTime()+bcTimePhaseCorr, res.GetChi2(), res.GetNdf() );
428 : }
429 58 : }// ALTRO
430 0 : else if ( fUseFALTRO )
431 : {// Fake ALTRO
432 0 : fTriggerRawDigitMaker->Add( bunchlist );
433 : }// Fake ALTRO
434 174 : } // End while over channel
435 : } // End while over DDL's, of input stream
436 :
437 4 : fTriggerRawDigitMaker->PostProcess();
438 :
439 4 : TrimDigits(digitsArr);
440 8 : }
441 :
442 : ///
443 : /// Remove entries with LGnoHG (unphysical), out of time window, and too bad chi2.
444 : ///
445 : void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr)
446 : {
447 : AliEMCALDigit *digit = 0;
448 :
449 : Int_t n = 0;
450 8 : Int_t nDigits = digitsArr->GetEntriesFast();
451 :
452 4 : TIter nextdigit(digitsArr);
453 :
454 126 : while ((digit = (AliEMCALDigit*) nextdigit()))
455 : {
456 57 : if (digit->GetType() == AliEMCALDigit::kLGnoHG)
457 : {
458 0 : AliDebug(1,Form("Remove digit with id %d, LGnoHG",digit->GetId()));
459 0 : digitsArr->Remove(digit);
460 : }
461 114 : else if(fTimeMin > digit->GetTime() || fTimeMax < digit->GetTime())
462 : {
463 0 : digitsArr->Remove(digit);
464 0 : AliDebug(1,Form("Remove digit with id %d, Bad Time %e",digit->GetId(), digit->GetTime()));
465 : }
466 57 : else if (0 > digit->GetChi2())
467 : {
468 0 : digitsArr->Remove(digit);
469 0 : AliDebug(1,Form("Remove digit with id %d, Bad Chi2 %e",digit->GetId(), digit->GetChi2()));
470 : }
471 : else
472 : {
473 57 : digit->SetIndexInList(n);
474 57 : n++;
475 : }
476 : }//while
477 :
478 4 : digitsArr->Compress();
479 :
480 20 : AliDebug(1,Form("N Digits before trimming : %d; after array compression %d",nDigits,digitsArr->GetEntriesFast()));
481 4 : }
482 :
483 : ///
484 : /// Select which fitting algo should be used.
485 : ///
486 : void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo)
487 : {
488 24 : delete fRawAnalyzer; // delete doesn't do anything if the pointer is 0x0
489 :
490 10 : fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer( fitAlgo );
491 10 : fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done, for the new methods
492 10 : fRawAnalyzer->SetOverflowCut ( OVERFLOWCUT );
493 10 : fRawAnalyzer->SetAmpCut(fNoiseThreshold);
494 10 : fRawAnalyzer->SetFitArrayCut(fNoiseThreshold);
495 10 : }
496 :
497 :
498 :
|
