adding ITS standalone macro (E.Biolcati)
[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) ;
379c5c09 247// time = CalibrateT(time,relId,lowGainFlag) ;
6f47f50d 248
249 //convert time from sample bin units to s
250 time*=fSampleToSec ;
379c5c09 251
252 if(energy <= 0.)
253 continue;
254
255 if (caloFlag == AliCaloRawStreamV3::kLowGain) {
256 new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
c2244d6d 257 if (sigLength>0 && fADCValuesLG!=0)
258 dynamic_cast<AliPHOSDigit*>(tmpLG.At(ilgDigit))->SetALTROSamplesLG(sigLength,fADCValuesLG);
379c5c09 259 ilgDigit++ ;
260 }
261 else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
262 if(fitter->IsOverflow()) //Keep this digit to replace it by Low Gain later.
263 //If there is no LogGain it wil be removed by cut on Min E
264 new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
265 else
266 new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
c2244d6d 267 if (sigLength>0 && fADCValuesHG!=0)
268 dynamic_cast<AliPHOSDigit*>(digits->At(iDigit))->SetALTROSamplesHG(sigLength,fADCValuesHG);
379c5c09 269 iDigit++;
270 }
271 } // End of NextChannel()
c0394bfb 272
379c5c09 273 //Now scan created LG and HG digits and keep only those which appeared in both lists
274 //replace energy of HighGain digits only if there is overflow
275 //negative energy (overflow)
276 digits->Sort() ;
277 tmpLG.Sort() ;
278 Int_t iLG = 0;
279 Int_t nLG1 = tmpLG.GetEntriesFast()-1 ;
280
281 for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) {
282 AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
283 if (!digHG) continue;
284 AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
285 while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
286 iLG++ ;
287 digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
288 }
289 absId=digHG->GetId() ;
290 fGeom->AbsToRelNumbering(absId,relId) ;
a7c985dc 291
379c5c09 292 if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
293 if(digHG->GetEnergy()<0.){ //This is overflow in HG
294 digHG->SetTime(digLG->GetTime()) ;
295 digHG->SetEnergy(digLG->GetEnergy()) ;
296 }
379c5c09 297 }
298 else{ //no pair - remove
fdd91c5a 299 if(digHG->GetEnergy()<0.) //no pair, in saturation
300 digits->RemoveAt(iDig) ;
88fb5e50 301 }
88fb5e50 302 }
379c5c09 303 } // End of NextDDL()
c0394bfb 304
305 CleanDigits(digits) ;
379c5c09 306
c0394bfb 307}
308//____________________________________________________________________________
309Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
310{
634d2178 311 // Convert EMC LG amplitude to HG (multipli by ~16)
312 // Calibration parameters are taken from calibration data base
c0394bfb 313 if(fgCalibData){
314 Int_t module = relId[0];
315 Int_t column = relId[3];
316 Int_t row = relId[2];
317 if(relId[1]==0) { // this is EMC
318 if(isLowGain){
319 amp*= fgCalibData->GetHighLowRatioEmc(module,column,row);
88fb5e50 320 }
c0394bfb 321 return amp ;
322 }
323 }
324 return 0;
325}
326//____________________________________________________________________________
327Double_t AliPHOSRawDigiProducer::CalibrateT(Double_t time, Int_t * relId, Bool_t /* isLowGain */)
328{
329 //Calibrate time
39569d36 330 time*=fPulseGenerator->GetRawFormatTimeTrigger() ;
c0394bfb 331 if(fgCalibData){
332 Int_t module = relId[0];
333 Int_t column = relId[3];
334 Int_t row = relId[2];
335 if(relId[1]==0) { // this is EMC
336 time += fgCalibData->GetTimeShiftEmc(module,column,row);
337 return time ;
338 }
339 }
340
341 return -999.;
342}
343//____________________________________________________________________________
344void AliPHOSRawDigiProducer::CleanDigits(TClonesArray * digits)
345{
346 // remove digits with amplitudes below threshold.
347 // remove digits in bad channels
348 // sort digits with icreasing AbsId
349
350 //remove digits in bad map and below threshold
351 Bool_t isBadMap = 0 ;
352 if(fgCalibData->GetNumOfEmcBadChannels()){
353 isBadMap=1 ;
354 }
355
356 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
357 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
358 if(!digit)
359 continue ;
360 if ( (IsInEMC(digit) && digit->GetEnergy() < fEmcMinE) ||
361 (IsInCPV(digit) && digit->GetEnergy() < fCpvMinE) ){
362 digits->RemoveAt(i) ;
363 continue ;
364 }
365 if(isBadMap){ //check bad map now
366 Int_t relid[4] ;
367 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
a65517aa 368 if(fgCalibData->IsBadChannelEmc(relid[0],relid[3],relid[2])){
c0394bfb 369 digits->RemoveAt(i) ;
88fb5e50 370 }
371 }
88fb5e50 372 }
373
c0394bfb 374 //Compress, sort and set indexes
375 digits->Compress() ;
376// digits->Sort(); already sorted earlier
377 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
378 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
379 digit->SetIndexInList(i) ;
380 }
381}
382//____________________________________________________________________________
383Bool_t AliPHOSRawDigiProducer::IsInEMC(AliPHOSDigit * digit) const
384{
385 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
386 return digit->GetId() <= fEmcCrystals ;
387
388}
389
390//____________________________________________________________________________
391Bool_t AliPHOSRawDigiProducer::IsInCPV(AliPHOSDigit * digit) const
392{
393 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
394 return digit->GetId() > fEmcCrystals ;
395}
396//____________________________________________________________________________
397void AliPHOSRawDigiProducer::GetCalibrationParameters()
398{
399 // Set calibration parameters:
400 // if calibration database exists, they are read from database,
401 // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
402 //
403 // It is a user responsilibity to open CDB before reconstruction, for example:
404 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
405
406 if (!fgCalibData){
407 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
408 }
409 if (fgCalibData->GetCalibDataEmc() == 0)
410 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
411 if (fgCalibData->GetCalibDataCpv() == 0)
412 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
88fb5e50 413}