Fix
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawDigiProducer.cxx
CommitLineData
7bc140a5 1/**************************************************************************
2 * Copyright(c) 2007, 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//This class produces PHOS digits of one event
379c5c09 19//using AliPHOSRawFitter.
7bc140a5 20//
21// For example:
22// TClonesArray *digits = new TClonesArray("AliPHOSDigit",100);
23// AliRawReader* rawReader = new AliRawReaderDate("2006run2211.raw");
24// AliPHOSRawDecoder dc(rawReader);
25// while (rawReader->NextEvent()) {
26// AliPHOSRawDigiProducer producer;
27// producer.MakeDigits(digits,&dc);
28// }
29
30// Author: Boris Polichtchouk
31
32// --- ROOT system ---
88fb5e50 33#include "TClonesArray.h"
34
7bc140a5 35// --- AliRoot header files ---
88fb5e50 36#include "AliPHOSRawDigiProducer.h"
379c5c09 37#include "AliPHOSRawFitterv0.h"
88fb5e50 38#include "AliPHOSGeometry.h"
39#include "AliPHOSDigit.h"
c0394bfb 40#include "AliPHOSCalibData.h"
39569d36 41#include "AliPHOSPulseGenerator.h"
379c5c09 42#include "AliCaloRawStreamV3.h"
c0394bfb 43#include "AliLog.h"
88fb5e50 44
45ClassImp(AliPHOSRawDigiProducer)
46
c0394bfb 47AliPHOSCalibData * AliPHOSRawDigiProducer::fgCalibData = 0 ;
48
49//--------------------------------------------------------------------------------------
7e88424f 50AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():
51 TObject(),
52 fEmcMinE(0.),
53 fCpvMinE(0.),
54 fSampleQualityCut(1.),
6f47f50d 55 fSampleToSec(0.),
7e88424f 56 fEmcCrystals(0),
57 fGeom(0),
379c5c09 58 fPulseGenerator(0),
59 fRawReader(0),
54ac223e 60 fRawStream(0),
61 fADCValuesLG(0),
62 fADCValuesHG(0)
7e88424f 63{
64 // Default constructor
70d93620 65
c0394bfb 66 fGeom=AliPHOSGeometry::GetInstance() ;
67 if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");
68
69 fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
39569d36 70 fPulseGenerator = new AliPHOSPulseGenerator();
c0394bfb 71 GetCalibrationParameters() ;
8be3a30a 72
c0394bfb 73}
379c5c09 74//--------------------------------------------------------------------------------------
75AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(AliRawReader *rawReader,
76 AliAltroMapping **mapping):
77 TObject(),
78 fEmcMinE(0.),
79 fCpvMinE(0.),
80 fSampleQualityCut(1.),
6f47f50d 81 fSampleToSec(0.),
379c5c09 82 fEmcCrystals(0),
83 fGeom(0),
84 fPulseGenerator(0),
85 fRawReader(rawReader),
54ac223e 86 fRawStream(0),
87 fADCValuesLG(0),
88 fADCValuesHG(0)
379c5c09 89{
90 // Default constructor
91
92 fGeom=AliPHOSGeometry::GetInstance() ;
93 if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");
94
95 fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
96 fPulseGenerator = new AliPHOSPulseGenerator();
97 GetCalibrationParameters() ;
98
99 fRawStream = new AliCaloRawStreamV3(rawReader,"PHOS",mapping);
100
101}
102//--------------------------------------------------------------------------------------
7e88424f 103AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):
104 TObject(),
105 fEmcMinE(0.),
106 fCpvMinE(0.),
107 fSampleQualityCut(1.),
6f47f50d 108 fSampleToSec(0.),
7e88424f 109 fEmcCrystals(0),
110 fGeom(0),
379c5c09 111 fPulseGenerator(0),
112 fRawReader(0),
54ac223e 113 fRawStream(0),
114 fADCValuesLG(0),
115 fADCValuesHG(0)
116
7e88424f 117{
118 // Copy constructor
c0394bfb 119
120 fEmcMinE = dp.fEmcMinE ;
121 fCpvMinE = dp.fCpvMinE ;
f78c9781 122 fSampleQualityCut = dp.fSampleQualityCut;
6f47f50d 123 fSampleToSec = dp.fSampleToSec ;
c0394bfb 124 fEmcCrystals = dp.fEmcCrystals ;
39569d36 125 fPulseGenerator = new AliPHOSPulseGenerator();
c0394bfb 126 fGeom = dp.fGeom ;
127}
128//--------------------------------------------------------------------------------------
7e88424f 129AliPHOSRawDigiProducer& AliPHOSRawDigiProducer::operator= (const AliPHOSRawDigiProducer &dp)
130{
131 // Assign operator
132
c0394bfb 133 if(&dp == this) return *this;
134
135 fEmcMinE = dp.fEmcMinE ;
136 fCpvMinE = dp.fCpvMinE ;
70d93620 137 fSampleQualityCut = dp.fSampleQualityCut ;
6f47f50d 138 fSampleToSec = dp.fSampleToSec ;
c0394bfb 139 fEmcCrystals = dp.fEmcCrystals ;
140 fGeom = dp.fGeom ;
39569d36 141 if(fPulseGenerator) delete fPulseGenerator ;
142 fPulseGenerator = new AliPHOSPulseGenerator();
c0394bfb 143 return *this;
144}
379c5c09 145//--------------------------------------------------------------------------------------
7e88424f 146AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer()
147{
379c5c09 148 // Destructor
7e88424f 149 if(fPulseGenerator) delete fPulseGenerator ;
150 fPulseGenerator=0 ;
379c5c09 151 delete fRawStream;
54ac223e 152 delete [] fADCValuesLG;
153 delete [] fADCValuesHG;
39569d36 154}
7bc140a5 155//--------------------------------------------------------------------------------------
379c5c09 156void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawFitterv0* fitter)
88fb5e50 157{
771123c7 158 // Create a temporary array of LG digits and then make digits from raw data
159
160 TClonesArray *tmpLG = new TClonesArray("AliPHOSDigit",10000) ;
161 MakeDigits(digits, tmpLG, fitter);
162 tmpLG->Delete();
163 delete tmpLG;
164}
165//--------------------------------------------------------------------------------------
166void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, TClonesArray *tmpDigLG, AliPHOSRawFitterv0* fitter)
167{
7bc140a5 168 //Makes the job.
771123c7 169 //TClonesArray *digits, *tmpDigLG and raw data fitter should be provided by calling function.
7bc140a5 170
88fb5e50 171 digits->Clear();
172
379c5c09 173 Int_t iDigit=0 ;
174 Int_t relId[4], absId=-1, caloFlag=-1;
fdd91c5a 175
c0394bfb 176 const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch.
6f47f50d 177
178 //Calculate conversion coeff. from Sample time step to seconds
179 //If OCDB contains negative or zero value - use one from RCU trailer
180 //Negative value in OCDB is used only for simulation of raw digits
181 if(fgCalibData->GetSampleTimeStep()>0.)
182 fSampleToSec=fgCalibData->GetSampleTimeStep() ;
183 else
184 fSampleToSec=fRawStream->GetTSample() ;
fdd91c5a 185
771123c7 186 // Clear a temporary array for LowGain digits
187 tmpDigLG->Clear();
c0394bfb 188 Int_t ilgDigit=0 ;
bafc1087 189
379c5c09 190 //Let fitter subtract pedestals in case of ZS
191 fitter->SetCalibData(fgCalibData) ;
c0394bfb 192
379c5c09 193 while (fRawStream->NextDDL()) {
194 while (fRawStream->NextChannel()) {
fc44a661 195 relId[0] = 5 - fRawStream->GetModule() ; // counts from 1 to 5
379c5c09 196 relId[1] = 0;
1281504a 197 relId[2] = fRawStream->GetCellX() + 1; // counts from 1 to 64
198 relId[3] = fRawStream->GetCellZ() + 1; // counts from 1 to 56
199 caloFlag = fRawStream->GetCaloFlag(); // 0=LG, 1=HG, 2=TRU
f4ca5622 200
201 if(caloFlag!=0 && caloFlag!=1) continue; //TRU data!
202
e10f780e 203 fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag);
fc44a661 204
c4dad9d2 205 if(fitter->GetAmpOffset()==0 && fitter->GetAmpThreshold()==0) {
206 short value = fRawStream->GetAltroCFG1();
207 bool ZeroSuppressionEnabled = (value >> 15) & 0x1;
208 if(ZeroSuppressionEnabled) {
209 short offset = (value >> 10) & 0xf;
210 short threshold = value & 0x3ff;
211 fitter->SubtractPedestals(kFALSE);
212 fitter->SetAmpOffset(offset);
213 fitter->SetAmpThreshold(threshold);
214 }
215 }
216
379c5c09 217 fGeom->RelToAbsNumbering(relId, absId);
f4ca5622 218
c01c5641 219 fitter->SetNBunches(0);
c2244d6d 220 Int_t sigStart =0 ;
221 Int_t sigLength=0 ;
222 while (fRawStream->NextBunch()) { //Take the first in time bunch
379c5c09 223 const UShort_t *sig = fRawStream->GetSignals();
54ac223e 224 sigStart = fRawStream->GetStartTimeBin();
225 sigLength = fRawStream->GetBunchLength();
92236b27 226 fitter->Eval(sig,sigStart,sigLength);
54ac223e 227 if (caloFlag == AliCaloRawStreamV3::kLowGain) {
228 delete [] fADCValuesLG;
229 fADCValuesLG = new Int_t[sigLength];
230 for (Int_t i=0; i<sigLength; i++)
231 fADCValuesLG[sigLength-i-1] = sig[i];
232 }
233 else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
234 delete [] fADCValuesHG;
235 fADCValuesHG = new Int_t[sigLength];
236 for (Int_t i=0; i<sigLength; i++)
237 fADCValuesHG[sigLength-i-1] = sig[i];
238 }
379c5c09 239 } // End of NextBunch()
70d93620 240
379c5c09 241
242 Double_t energy = fitter->GetEnergy() ;
243 Double_t time = fitter->GetTime() ;
244 if(energy<=baseLine) //in ADC channels
245 continue ;
88fb5e50 246
379c5c09 247 //remove digits with bad shape. Fitter should calculate quality so that
248 //in default case quality [0,1], while larger values of quality mean somehow
249 //corrupted samples, 999 means obviously corrupted sample.
250 //It is difficult to fit samples with overflow (even setting cut on overflow values)
251 //because too few points are left to fit. So we do not evaluate samples with overflow
c0394bfb 252
379c5c09 253 if(fitter->GetSignalQuality() > fSampleQualityCut && !(fitter->IsOverflow()))
254 continue ;
255
817744df 256 energy = CalibrateE(energy,relId,!caloFlag) ;
6f47f50d 257
258 //convert time from sample bin units to s
259 time*=fSampleToSec ;
b879fc0b 260//CalibrateT moved to Clusterizer
261// time = CalibrateT(time,relId,!caloFlag) ;
262 // subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0:
911980a0 263 time -= fRawStream->GetL1Phase();
b879fc0b 264
265
379c5c09 266
267 if(energy <= 0.)
268 continue;
269
270 if (caloFlag == AliCaloRawStreamV3::kLowGain) {
771123c7 271 new((*tmpDigLG)[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
c2244d6d 272 if (sigLength>0 && fADCValuesLG!=0)
09828f52 273 static_cast<AliPHOSDigit*>(tmpDigLG->At(ilgDigit))->SetALTROSamplesLG(sigLength,fADCValuesLG);
379c5c09 274 ilgDigit++ ;
275 }
276 else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
277 if(fitter->IsOverflow()) //Keep this digit to replace it by Low Gain later.
278 //If there is no LogGain it wil be removed by cut on Min E
279 new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
280 else
281 new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
c2244d6d 282 if (sigLength>0 && fADCValuesHG!=0)
09828f52 283 static_cast<AliPHOSDigit*>(digits->At(iDigit))->SetALTROSamplesHG(sigLength,fADCValuesHG);
379c5c09 284 iDigit++;
285 }
286 } // End of NextChannel()
c0394bfb 287
379c5c09 288 //Now scan created LG and HG digits and keep only those which appeared in both lists
289 //replace energy of HighGain digits only if there is overflow
290 //negative energy (overflow)
291 digits->Sort() ;
771123c7 292 tmpDigLG->Sort() ;
379c5c09 293 Int_t iLG = 0;
771123c7 294 Int_t nLG1 = tmpDigLG->GetEntriesFast()-1 ;
379c5c09 295
296 for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) {
297 AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
298 if (!digHG) continue;
771123c7 299 AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpDigLG->At(iLG)) ;
379c5c09 300 while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
301 iLG++ ;
771123c7 302 digLG = dynamic_cast<AliPHOSDigit*>(tmpDigLG->At(iLG)) ;
379c5c09 303 }
304 absId=digHG->GetId() ;
305 fGeom->AbsToRelNumbering(absId,relId) ;
a7c985dc 306
379c5c09 307 if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
308 if(digHG->GetEnergy()<0.){ //This is overflow in HG
309 digHG->SetTime(digLG->GetTime()) ;
310 digHG->SetEnergy(digLG->GetEnergy()) ;
311 }
379c5c09 312 }
313 else{ //no pair - remove
fdd91c5a 314 if(digHG->GetEnergy()<0.) //no pair, in saturation
315 digits->RemoveAt(iDig) ;
88fb5e50 316 }
88fb5e50 317 }
379c5c09 318 } // End of NextDDL()
c0394bfb 319
320 CleanDigits(digits) ;
379c5c09 321
c0394bfb 322}
323//____________________________________________________________________________
324Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
325{
634d2178 326 // Convert EMC LG amplitude to HG (multipli by ~16)
327 // Calibration parameters are taken from calibration data base
c0394bfb 328 if(fgCalibData){
329 Int_t module = relId[0];
330 Int_t column = relId[3];
331 Int_t row = relId[2];
332 if(relId[1]==0) { // this is EMC
333 if(isLowGain){
334 amp*= fgCalibData->GetHighLowRatioEmc(module,column,row);
88fb5e50 335 }
c0394bfb 336 return amp ;
337 }
338 }
339 return 0;
340}
341//____________________________________________________________________________
342Double_t AliPHOSRawDigiProducer::CalibrateT(Double_t time, Int_t * relId, Bool_t /* isLowGain */)
343{
344 //Calibrate time
345 if(fgCalibData){
346 Int_t module = relId[0];
347 Int_t column = relId[3];
348 Int_t row = relId[2];
349 if(relId[1]==0) { // this is EMC
350 time += fgCalibData->GetTimeShiftEmc(module,column,row);
351 return time ;
352 }
353 }
354
355 return -999.;
356}
357//____________________________________________________________________________
358void AliPHOSRawDigiProducer::CleanDigits(TClonesArray * digits)
359{
360 // remove digits with amplitudes below threshold.
361 // remove digits in bad channels
362 // sort digits with icreasing AbsId
363
364 //remove digits in bad map and below threshold
365 Bool_t isBadMap = 0 ;
366 if(fgCalibData->GetNumOfEmcBadChannels()){
367 isBadMap=1 ;
368 }
369
370 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
371 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
372 if(!digit)
373 continue ;
374 if ( (IsInEMC(digit) && digit->GetEnergy() < fEmcMinE) ||
375 (IsInCPV(digit) && digit->GetEnergy() < fCpvMinE) ){
376 digits->RemoveAt(i) ;
377 continue ;
378 }
379 if(isBadMap){ //check bad map now
380 Int_t relid[4] ;
381 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
a65517aa 382 if(fgCalibData->IsBadChannelEmc(relid[0],relid[3],relid[2])){
c0394bfb 383 digits->RemoveAt(i) ;
88fb5e50 384 }
385 }
88fb5e50 386 }
387
c0394bfb 388 //Compress, sort and set indexes
389 digits->Compress() ;
390// digits->Sort(); already sorted earlier
391 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
392 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
393 digit->SetIndexInList(i) ;
394 }
395}
396//____________________________________________________________________________
397Bool_t AliPHOSRawDigiProducer::IsInEMC(AliPHOSDigit * digit) const
398{
399 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
400 return digit->GetId() <= fEmcCrystals ;
401
402}
403
404//____________________________________________________________________________
405Bool_t AliPHOSRawDigiProducer::IsInCPV(AliPHOSDigit * digit) const
406{
407 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
408 return digit->GetId() > fEmcCrystals ;
409}
410//____________________________________________________________________________
411void AliPHOSRawDigiProducer::GetCalibrationParameters()
412{
413 // Set calibration parameters:
414 // if calibration database exists, they are read from database,
415 // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
416 //
417 // It is a user responsilibity to open CDB before reconstruction, for example:
418 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
419
420 if (!fgCalibData){
421 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
422 }
423 if (fgCalibData->GetCalibDataEmc() == 0)
424 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
425 if (fgCalibData->GetCalibDataCpv() == 0)
426 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
88fb5e50 427}