andrew updates from sat night.
[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{
7bc140a5 158 //Makes the job.
379c5c09 159 //TClonesArray *digits and raw data fitter should be provided by calling function.
7bc140a5 160
88fb5e50 161 digits->Clear();
162
379c5c09 163 Int_t iDigit=0 ;
164 Int_t relId[4], absId=-1, caloFlag=-1;
fdd91c5a 165
c0394bfb 166 const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch.
6f47f50d 167
168 //Calculate conversion coeff. from Sample time step to seconds
169 //If OCDB contains negative or zero value - use one from RCU trailer
170 //Negative value in OCDB is used only for simulation of raw digits
171 if(fgCalibData->GetSampleTimeStep()>0.)
172 fSampleToSec=fgCalibData->GetSampleTimeStep() ;
173 else
174 fSampleToSec=fRawStream->GetTSample() ;
fdd91c5a 175
c0394bfb 176 //Temporary array for LowGain digits
177 TClonesArray tmpLG("AliPHOSDigit",10000) ;
178 Int_t ilgDigit=0 ;
bafc1087 179
379c5c09 180 //Let fitter subtract pedestals in case of ZS
181 fitter->SetCalibData(fgCalibData) ;
c0394bfb 182
379c5c09 183 while (fRawStream->NextDDL()) {
184 while (fRawStream->NextChannel()) {
fc44a661 185 relId[0] = 5 - fRawStream->GetModule() ; // counts from 1 to 5
379c5c09 186 relId[1] = 0;
1281504a 187 relId[2] = fRawStream->GetCellX() + 1; // counts from 1 to 64
188 relId[3] = fRawStream->GetCellZ() + 1; // counts from 1 to 56
189 caloFlag = fRawStream->GetCaloFlag(); // 0=LG, 1=HG, 2=TRU
f4ca5622 190
191 if(caloFlag!=0 && caloFlag!=1) continue; //TRU data!
192
e10f780e 193 fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag);
fc44a661 194
c4dad9d2 195 if(fitter->GetAmpOffset()==0 && fitter->GetAmpThreshold()==0) {
196 short value = fRawStream->GetAltroCFG1();
197 bool ZeroSuppressionEnabled = (value >> 15) & 0x1;
198 if(ZeroSuppressionEnabled) {
199 short offset = (value >> 10) & 0xf;
200 short threshold = value & 0x3ff;
201 fitter->SubtractPedestals(kFALSE);
202 fitter->SetAmpOffset(offset);
203 fitter->SetAmpThreshold(threshold);
204 }
205 }
206
379c5c09 207 fGeom->RelToAbsNumbering(relId, absId);
f4ca5622 208
c01c5641 209 fitter->SetNBunches(0);
c2244d6d 210 Int_t sigStart =0 ;
211 Int_t sigLength=0 ;
212 while (fRawStream->NextBunch()) { //Take the first in time bunch
379c5c09 213 const UShort_t *sig = fRawStream->GetSignals();
54ac223e 214 sigStart = fRawStream->GetStartTimeBin();
215 sigLength = fRawStream->GetBunchLength();
92236b27 216 fitter->Eval(sig,sigStart,sigLength);
54ac223e 217 if (caloFlag == AliCaloRawStreamV3::kLowGain) {
218 delete [] fADCValuesLG;
219 fADCValuesLG = new Int_t[sigLength];
220 for (Int_t i=0; i<sigLength; i++)
221 fADCValuesLG[sigLength-i-1] = sig[i];
222 }
223 else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
224 delete [] fADCValuesHG;
225 fADCValuesHG = new Int_t[sigLength];
226 for (Int_t i=0; i<sigLength; i++)
227 fADCValuesHG[sigLength-i-1] = sig[i];
228 }
379c5c09 229 } // End of NextBunch()
70d93620 230
379c5c09 231
232 Double_t energy = fitter->GetEnergy() ;
233 Double_t time = fitter->GetTime() ;
234 if(energy<=baseLine) //in ADC channels
235 continue ;
88fb5e50 236
379c5c09 237 //remove digits with bad shape. Fitter should calculate quality so that
238 //in default case quality [0,1], while larger values of quality mean somehow
239 //corrupted samples, 999 means obviously corrupted sample.
240 //It is difficult to fit samples with overflow (even setting cut on overflow values)
241 //because too few points are left to fit. So we do not evaluate samples with overflow
c0394bfb 242
379c5c09 243 if(fitter->GetSignalQuality() > fSampleQualityCut && !(fitter->IsOverflow()))
244 continue ;
245
817744df 246 energy = CalibrateE(energy,relId,!caloFlag) ;
6f47f50d 247
248 //convert time from sample bin units to s
249 time*=fSampleToSec ;
b879fc0b 250//CalibrateT moved to Clusterizer
251// time = CalibrateT(time,relId,!caloFlag) ;
252 // subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0:
95476bb5 253 time += fRawStream->GetL1Phase();
b879fc0b 254
255
379c5c09 256
257 if(energy <= 0.)
258 continue;
259
260 if (caloFlag == AliCaloRawStreamV3::kLowGain) {
261 new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
c2244d6d 262 if (sigLength>0 && fADCValuesLG!=0)
263 dynamic_cast<AliPHOSDigit*>(tmpLG.At(ilgDigit))->SetALTROSamplesLG(sigLength,fADCValuesLG);
379c5c09 264 ilgDigit++ ;
265 }
266 else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
267 if(fitter->IsOverflow()) //Keep this digit to replace it by Low Gain later.
268 //If there is no LogGain it wil be removed by cut on Min E
269 new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
270 else
271 new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
c2244d6d 272 if (sigLength>0 && fADCValuesHG!=0)
273 dynamic_cast<AliPHOSDigit*>(digits->At(iDigit))->SetALTROSamplesHG(sigLength,fADCValuesHG);
379c5c09 274 iDigit++;
275 }
276 } // End of NextChannel()
c0394bfb 277
379c5c09 278 //Now scan created LG and HG digits and keep only those which appeared in both lists
279 //replace energy of HighGain digits only if there is overflow
280 //negative energy (overflow)
281 digits->Sort() ;
282 tmpLG.Sort() ;
283 Int_t iLG = 0;
284 Int_t nLG1 = tmpLG.GetEntriesFast()-1 ;
285
286 for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) {
287 AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
288 if (!digHG) continue;
289 AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
290 while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
291 iLG++ ;
292 digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
293 }
294 absId=digHG->GetId() ;
295 fGeom->AbsToRelNumbering(absId,relId) ;
a7c985dc 296
379c5c09 297 if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
298 if(digHG->GetEnergy()<0.){ //This is overflow in HG
299 digHG->SetTime(digLG->GetTime()) ;
300 digHG->SetEnergy(digLG->GetEnergy()) ;
301 }
379c5c09 302 }
303 else{ //no pair - remove
fdd91c5a 304 if(digHG->GetEnergy()<0.) //no pair, in saturation
305 digits->RemoveAt(iDig) ;
88fb5e50 306 }
88fb5e50 307 }
379c5c09 308 } // End of NextDDL()
c0394bfb 309
310 CleanDigits(digits) ;
379c5c09 311
c0394bfb 312}
313//____________________________________________________________________________
314Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
315{
634d2178 316 // Convert EMC LG amplitude to HG (multipli by ~16)
317 // Calibration parameters are taken from calibration data base
c0394bfb 318 if(fgCalibData){
319 Int_t module = relId[0];
320 Int_t column = relId[3];
321 Int_t row = relId[2];
322 if(relId[1]==0) { // this is EMC
323 if(isLowGain){
324 amp*= fgCalibData->GetHighLowRatioEmc(module,column,row);
88fb5e50 325 }
c0394bfb 326 return amp ;
327 }
328 }
329 return 0;
330}
331//____________________________________________________________________________
332Double_t AliPHOSRawDigiProducer::CalibrateT(Double_t time, Int_t * relId, Bool_t /* isLowGain */)
333{
334 //Calibrate time
335 if(fgCalibData){
336 Int_t module = relId[0];
337 Int_t column = relId[3];
338 Int_t row = relId[2];
339 if(relId[1]==0) { // this is EMC
340 time += fgCalibData->GetTimeShiftEmc(module,column,row);
341 return time ;
342 }
343 }
344
345 return -999.;
346}
347//____________________________________________________________________________
348void AliPHOSRawDigiProducer::CleanDigits(TClonesArray * digits)
349{
350 // remove digits with amplitudes below threshold.
351 // remove digits in bad channels
352 // sort digits with icreasing AbsId
353
354 //remove digits in bad map and below threshold
355 Bool_t isBadMap = 0 ;
356 if(fgCalibData->GetNumOfEmcBadChannels()){
357 isBadMap=1 ;
358 }
359
360 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
361 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
362 if(!digit)
363 continue ;
364 if ( (IsInEMC(digit) && digit->GetEnergy() < fEmcMinE) ||
365 (IsInCPV(digit) && digit->GetEnergy() < fCpvMinE) ){
366 digits->RemoveAt(i) ;
367 continue ;
368 }
369 if(isBadMap){ //check bad map now
370 Int_t relid[4] ;
371 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
a65517aa 372 if(fgCalibData->IsBadChannelEmc(relid[0],relid[3],relid[2])){
c0394bfb 373 digits->RemoveAt(i) ;
88fb5e50 374 }
375 }
88fb5e50 376 }
377
c0394bfb 378 //Compress, sort and set indexes
379 digits->Compress() ;
380// digits->Sort(); already sorted earlier
381 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
382 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
383 digit->SetIndexInList(i) ;
384 }
385}
386//____________________________________________________________________________
387Bool_t AliPHOSRawDigiProducer::IsInEMC(AliPHOSDigit * digit) const
388{
389 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
390 return digit->GetId() <= fEmcCrystals ;
391
392}
393
394//____________________________________________________________________________
395Bool_t AliPHOSRawDigiProducer::IsInCPV(AliPHOSDigit * digit) const
396{
397 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
398 return digit->GetId() > fEmcCrystals ;
399}
400//____________________________________________________________________________
401void AliPHOSRawDigiProducer::GetCalibrationParameters()
402{
403 // Set calibration parameters:
404 // if calibration database exists, they are read from database,
405 // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
406 //
407 // It is a user responsilibity to open CDB before reconstruction, for example:
408 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
409
410 if (!fgCalibData){
411 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
412 }
413 if (fgCalibData->GetCalibDataEmc() == 0)
414 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
415 if (fgCalibData->GetCalibDataCpv() == 0)
416 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
88fb5e50 417}