[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawDigiProducer.cxx

/**************************************************************************
 * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//This class produces PHOS digits of one event
//using AliPHOSRawFitter. 
//
//   For example:
//   TClonesArray *digits = new TClonesArray("AliPHOSDigit",100);
//   AliRawReader* rawReader = new AliRawReaderDate("2006run2211.raw");
//   AliPHOSRawDecoder dc(rawReader);
//   while (rawReader->NextEvent()) {
//     AliPHOSRawDigiProducer producer;
//     producer.MakeDigits(digits,&dc);
//   }

// Author: Boris Polichtchouk

// --- ROOT system ---
#include "TClonesArray.h"

// --- AliRoot header files ---
#include "AliPHOSRawDigiProducer.h"
#include "AliPHOSRawFitterv0.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSDigit.h"
#include "AliPHOSCalibData.h"
#include "AliPHOSPulseGenerator.h"
#include "AliCaloRawStreamV3.h"
#include "AliLog.h"

ClassImp(AliPHOSRawDigiProducer)

AliPHOSCalibData * AliPHOSRawDigiProducer::fgCalibData  = 0 ; 

//--------------------------------------------------------------------------------------
AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():
  TObject(),
  fEmcMinE(0.),
  fCpvMinE(0.),
  fSampleQualityCut(1.),
  fSampleToSec(0.),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0),
  fRawReader(0),
  fRawStream(0),
  fADCValuesLG(0),
  fADCValuesHG(0)
{
  // Default constructor

  fGeom=AliPHOSGeometry::GetInstance() ;
  if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");

  fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
  fPulseGenerator = new AliPHOSPulseGenerator();
  GetCalibrationParameters() ; 

}
//--------------------------------------------------------------------------------------
AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(AliRawReader *rawReader,
					       AliAltroMapping **mapping):
  TObject(),
  fEmcMinE(0.),
  fCpvMinE(0.),
  fSampleQualityCut(1.),
  fSampleToSec(0.),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0),
  fRawReader(rawReader),
  fRawStream(0),
  fADCValuesLG(0),
  fADCValuesHG(0)
{
  // Default constructor

  fGeom=AliPHOSGeometry::GetInstance() ;
  if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");

  fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
  fPulseGenerator = new AliPHOSPulseGenerator();
  GetCalibrationParameters() ; 

  fRawStream = new AliCaloRawStreamV3(rawReader,"PHOS",mapping);

}
//--------------------------------------------------------------------------------------
AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):
  TObject(),
  fEmcMinE(0.),
  fCpvMinE(0.),
  fSampleQualityCut(1.),
  fSampleToSec(0.),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0),
  fRawReader(0),
  fRawStream(0),
  fADCValuesLG(0),
  fADCValuesHG(0)

{                                                          
  // Copy constructor

  fEmcMinE = dp.fEmcMinE ;
  fCpvMinE = dp.fCpvMinE ;
  fSampleQualityCut = dp.fSampleQualityCut;
  fSampleToSec = dp.fSampleToSec ;
  fEmcCrystals = dp.fEmcCrystals ;
  fPulseGenerator = new AliPHOSPulseGenerator();
  fGeom = dp.fGeom ;
}
//--------------------------------------------------------------------------------------
AliPHOSRawDigiProducer& AliPHOSRawDigiProducer::operator= (const AliPHOSRawDigiProducer &dp)
{
  // Assign operator

  if(&dp == this) return *this;

  fEmcMinE = dp.fEmcMinE ;
  fCpvMinE = dp.fCpvMinE ;
  fSampleQualityCut = dp.fSampleQualityCut ;
  fSampleToSec = dp.fSampleToSec ;
  fEmcCrystals = dp.fEmcCrystals ;
  fGeom = dp.fGeom ;
  if(fPulseGenerator) delete fPulseGenerator ;
  fPulseGenerator = new AliPHOSPulseGenerator();
  return  *this;
} 
//--------------------------------------------------------------------------------------
AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer()
{
  // Destructor
  if(fPulseGenerator) delete fPulseGenerator ;
  fPulseGenerator=0 ;
  delete fRawStream;
  delete [] fADCValuesLG;
  delete [] fADCValuesHG;
}
//--------------------------------------------------------------------------------------
void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawFitterv0* fitter) 
{
  // Create a temporary array of LG digits and then make digits from raw data

  TClonesArray *tmpLG = new TClonesArray("AliPHOSDigit",10000) ;
  MakeDigits(digits, tmpLG, fitter);
  tmpLG->Delete();
  delete tmpLG;
}
//--------------------------------------------------------------------------------------
void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, TClonesArray *tmpDigLG, AliPHOSRawFitterv0* fitter) 
{
  //Makes the job.
  //TClonesArray *digits, *tmpDigLG and raw data fitter should be provided by calling function.

  digits->Clear();
 
  Int_t iDigit=0 ;
  Int_t relId[4], absId=-1, caloFlag=-1;
  
  const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch. 

  //Calculate conversion coeff. from Sample time step to seconds
  //If OCDB contains negative or zero value - use one from RCU trailer
  //Negative value in OCDB is used only for simulation of raw digits
  if(fgCalibData->GetSampleTimeStep()>0.)
    fSampleToSec=fgCalibData->GetSampleTimeStep() ;
  else
    fSampleToSec=fRawStream->GetTSample() ;
  
  // Clear a temporary array for LowGain digits
  tmpDigLG->Clear();
  Int_t ilgDigit=0 ;

  //Let fitter subtract pedestals in case of ZS
  fitter->SetCalibData(fgCalibData) ;
  
  while (fRawStream->NextDDL()) {
    while (fRawStream->NextChannel()) {
      relId[0] = 5 - fRawStream->GetModule() ; // counts from 1 to 5
      relId[1] = 0;
      relId[2] = fRawStream->GetCellX()  + 1; // counts from 1 to 64
      relId[3] = fRawStream->GetCellZ()  + 1; // counts from 1 to 56
      caloFlag = fRawStream->GetCaloFlag();   // 0=LG, 1=HG, 2=TRU
      
      if(caloFlag!=0 && caloFlag!=1) continue; //TRU data!
      
      fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag);

      if(fitter->GetAmpOffset()==0 && fitter->GetAmpThreshold()==0) {
	short value = fRawStream->GetAltroCFG1();
	bool ZeroSuppressionEnabled = (value >> 15) & 0x1;
	if(ZeroSuppressionEnabled) {
	  short offset = (value >> 10) & 0xf;
	  short threshold = value & 0x3ff;
	  fitter->SubtractPedestals(kFALSE);
	  fitter->SetAmpOffset(offset);
	  fitter->SetAmpThreshold(threshold);
	}
      }
      
      fGeom->RelToAbsNumbering(relId, absId);
      
      fitter->SetNBunches(0);
      Int_t sigStart =0 ;
      Int_t sigLength=0 ;
      while (fRawStream->NextBunch()) { //Take the first in time bunch
	const UShort_t *sig = fRawStream->GetSignals();
	sigStart  = fRawStream->GetStartTimeBin();
	sigLength = fRawStream->GetBunchLength();
	fitter->Eval(sig,sigStart,sigLength);
	if      (caloFlag == AliCaloRawStreamV3::kLowGain) {
	  delete [] fADCValuesLG;
	  fADCValuesLG = new Int_t[sigLength];
	  for (Int_t i=0; i<sigLength; i++)
	    fADCValuesLG[sigLength-i-1] = sig[i];
	}
	else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
	  delete [] fADCValuesHG;
	  fADCValuesHG = new Int_t[sigLength];
	  for (Int_t i=0; i<sigLength; i++)
	    fADCValuesHG[sigLength-i-1] = sig[i];
	}
      } // End of NextBunch()

      
      Double_t energy = fitter->GetEnergy() ; 
      Double_t time   = fitter->GetTime() ;
      if(energy<=baseLine) //in ADC channels
	continue ;

      //remove digits with bad shape. Fitter should calculate quality so that 
      //in default case quality [0,1], while larger values of quality mean somehow 
      //corrupted samples, 999 means obviously corrupted sample.
      //It is difficult to fit samples with overflow (even setting cut on overflow values)
      //because too few points are left to fit. So we do not evaluate samples with overflow

      if(fitter->GetSignalQuality() > fSampleQualityCut && !(fitter->IsOverflow()))
	continue ;
      
      energy = CalibrateE(energy,relId,!caloFlag) ;

      //convert time from sample bin units to s
      time*=fSampleToSec ;
//CalibrateT moved to Clusterizer
//      time = CalibrateT(time,relId,!caloFlag) ;
      // subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0:
      time -= fRawStream->GetL1Phase();


      if(energy <= 0.) 
	continue;
      
      if (caloFlag == AliCaloRawStreamV3::kLowGain) {
	new((*tmpDigLG)[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
	if (sigLength>0 && fADCValuesLG!=0)
	  static_cast<AliPHOSDigit*>(tmpDigLG->At(ilgDigit))->SetALTROSamplesLG(sigLength,fADCValuesLG);
	ilgDigit++ ; 
      }
      else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
	if(fitter->IsOverflow()) //Keep this digit to replace it by Low Gain later.
	  //If there is no LogGain it wil be removed by cut on Min E
	  new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
	else
	  new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
	if (sigLength>0 && fADCValuesHG!=0)
	  static_cast<AliPHOSDigit*>(digits->At(iDigit))->SetALTROSamplesHG(sigLength,fADCValuesHG);
	iDigit++;
      }
    } // End of NextChannel()

    //Now scan created LG and HG digits and keep only those which appeared in both lists 
    //replace energy of HighGain digits only if there is overflow
    //negative energy (overflow)
    digits->Sort() ;
    tmpDigLG->Sort() ;
    Int_t iLG = 0;
    Int_t nLG1 = tmpDigLG->GetEntriesFast()-1 ;
    
    for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) { 
      AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
      if (!digHG) continue;
      AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpDigLG->At(iLG)) ;
      while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
	iLG++ ;
	digLG = dynamic_cast<AliPHOSDigit*>(tmpDigLG->At(iLG)) ;
      }
      absId=digHG->GetId() ;
      fGeom->AbsToRelNumbering(absId,relId) ;
 
      if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
	if(digHG->GetEnergy()<0.){ //This is overflow in HG
	  digHG->SetTime(digLG->GetTime()) ;
	  digHG->SetEnergy(digLG->GetEnergy()) ;
	}
      }
      else{ //no pair - remove
	if(digHG->GetEnergy()<0.) //no pair, in saturation
	  digits->RemoveAt(iDig) ;                                                          
      }
    }
  } // End of NextDDL()

  CleanDigits(digits) ;
  
}
//____________________________________________________________________________
Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
{
  // Convert EMC LG amplitude to HG (multipli by ~16)
  // Calibration parameters are taken from calibration data base 
  if(fgCalibData){ 
    Int_t   module = relId[0];  
    Int_t   column = relId[3];
    Int_t   row    = relId[2];
    if(relId[1]==0) { // this is EMC 
      if(isLowGain){
        amp*= fgCalibData->GetHighLowRatioEmc(module,column,row);
      }
      return amp ;         
    }         
  }          
  return 0;        
}
//____________________________________________________________________________
Double_t AliPHOSRawDigiProducer::CalibrateT(Double_t time, Int_t * relId, Bool_t /* isLowGain */)
{
  //Calibrate time
  if(fgCalibData){
    Int_t   module = relId[0];
    Int_t   column = relId[3];
    Int_t   row    = relId[2];
    if(relId[1]==0) { // this is EMC
      time += fgCalibData->GetTimeShiftEmc(module,column,row);                   
      return time ;             
    }
  }
 
  return -999.;
}
//____________________________________________________________________________
void AliPHOSRawDigiProducer::CleanDigits(TClonesArray * digits)
{
  // remove digits with amplitudes below threshold.
  // remove digits in bad channels
  // sort digits with icreasing AbsId
  
  //remove digits in bad map and below threshold
  Bool_t isBadMap = 0 ;
  if(fgCalibData->GetNumOfEmcBadChannels()){
    isBadMap=1 ;
  }
  
  for(Int_t i=0; i<digits->GetEntriesFast(); i++){
    AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
    if(!digit)
      continue  ;
    if ( (IsInEMC(digit) && digit->GetEnergy() < fEmcMinE) ||
	 (IsInCPV(digit) && digit->GetEnergy() < fCpvMinE) ){
      digits->RemoveAt(i) ;
      continue ;
    }
    if(isBadMap){ //check bad map now
      Int_t relid[4] ;
      fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
      if(fgCalibData->IsBadChannelEmc(relid[0],relid[3],relid[2])){
	digits->RemoveAt(i) ;
      }
    }
  }

  //Compress, sort and set indexes
  digits->Compress() ;
//  digits->Sort(); already sorted earlier
  for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) { 
    AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ; 
    digit->SetIndexInList(i) ;     
  }
}
//____________________________________________________________________________
Bool_t AliPHOSRawDigiProducer::IsInEMC(AliPHOSDigit * digit) const
{
  // Tells if (true) or not (false) the digit is in a PHOS-EMC module
  return digit->GetId() <= fEmcCrystals ;

}

//____________________________________________________________________________
Bool_t AliPHOSRawDigiProducer::IsInCPV(AliPHOSDigit * digit) const
{
  // Tells if (true) or not (false) the digit is in a PHOS-CPV module
  return digit->GetId() > fEmcCrystals ;
}
//____________________________________________________________________________
void AliPHOSRawDigiProducer::GetCalibrationParameters() 
{
  // Set calibration parameters:
  // if calibration database exists, they are read from database,
  // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
  //
  // It is a user responsilibity to open CDB before reconstruction, for example: 
  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");

  if (!fgCalibData){
    fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
  }
  if (fgCalibData->GetCalibDataEmc() == 0)
    AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
  if (fgCalibData->GetCalibDataCpv() == 0)
    AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
}
Commit	Line	Data
	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
	19	//using AliPHOSRawFitter.
	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 ---
	33	#include "TClonesArray.h"
	34
	35	// --- AliRoot header files ---
	36	#include "AliPHOSRawDigiProducer.h"
	37	#include "AliPHOSRawFitterv0.h"
	38	#include "AliPHOSGeometry.h"
	39	#include "AliPHOSDigit.h"
	40	#include "AliPHOSCalibData.h"
	41	#include "AliPHOSPulseGenerator.h"
	42	#include "AliCaloRawStreamV3.h"
	43	#include "AliLog.h"
	44
	45	ClassImp(AliPHOSRawDigiProducer)
	46
	47	AliPHOSCalibData * AliPHOSRawDigiProducer::fgCalibData = 0 ;
	48
	49	//--------------------------------------------------------------------------------------
	50	AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():
	51	TObject(),
	52	fEmcMinE(0.),
	53	fCpvMinE(0.),
	54	fSampleQualityCut(1.),
	55	fSampleToSec(0.),
	56	fEmcCrystals(0),
	57	fGeom(0),
	58	fPulseGenerator(0),
	59	fRawReader(0),
	60	fRawStream(0),
	61	fADCValuesLG(0),
	62	fADCValuesHG(0)
	63	{
	64	// Default constructor
	65
	66	fGeom=AliPHOSGeometry::GetInstance() ;
	67	if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");
	68
	69	fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
	70	fPulseGenerator = new AliPHOSPulseGenerator();
	71	GetCalibrationParameters() ;
	72
	73	}
	74	//--------------------------------------------------------------------------------------
	75	AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(AliRawReader *rawReader,
	76	AliAltroMapping **mapping):
	77	TObject(),
	78	fEmcMinE(0.),
	79	fCpvMinE(0.),
	80	fSampleQualityCut(1.),
	81	fSampleToSec(0.),
	82	fEmcCrystals(0),
	83	fGeom(0),
	84	fPulseGenerator(0),
	85	fRawReader(rawReader),
	86	fRawStream(0),
	87	fADCValuesLG(0),
	88	fADCValuesHG(0)
	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	//--------------------------------------------------------------------------------------
	103	AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):
	104	TObject(),
	105	fEmcMinE(0.),
	106	fCpvMinE(0.),
	107	fSampleQualityCut(1.),
	108	fSampleToSec(0.),
	109	fEmcCrystals(0),
	110	fGeom(0),
	111	fPulseGenerator(0),
	112	fRawReader(0),
	113	fRawStream(0),
	114	fADCValuesLG(0),
	115	fADCValuesHG(0)
	116
	117	{
	118	// Copy constructor
	119
	120	fEmcMinE = dp.fEmcMinE ;
	121	fCpvMinE = dp.fCpvMinE ;
	122	fSampleQualityCut = dp.fSampleQualityCut;
	123	fSampleToSec = dp.fSampleToSec ;
	124	fEmcCrystals = dp.fEmcCrystals ;
	125	fPulseGenerator = new AliPHOSPulseGenerator();
	126	fGeom = dp.fGeom ;
	127	}
	128	//--------------------------------------------------------------------------------------
	129	AliPHOSRawDigiProducer& AliPHOSRawDigiProducer::operator= (const AliPHOSRawDigiProducer &dp)
	130	{
	131	// Assign operator
	132
	133	if(&dp == this) return *this;
	134
	135	fEmcMinE = dp.fEmcMinE ;
	136	fCpvMinE = dp.fCpvMinE ;
	137	fSampleQualityCut = dp.fSampleQualityCut ;
	138	fSampleToSec = dp.fSampleToSec ;
	139	fEmcCrystals = dp.fEmcCrystals ;
	140	fGeom = dp.fGeom ;
	141	if(fPulseGenerator) delete fPulseGenerator ;
	142	fPulseGenerator = new AliPHOSPulseGenerator();
	143	return *this;
	144	}
	145	//--------------------------------------------------------------------------------------
	146	AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer()
	147	{
	148	// Destructor
	149	if(fPulseGenerator) delete fPulseGenerator ;
	150	fPulseGenerator=0 ;
	151	delete fRawStream;
	152	delete [] fADCValuesLG;
	153	delete [] fADCValuesHG;
	154	}
	155	//--------------------------------------------------------------------------------------
	156	void AliPHOSRawDigiProducer::MakeDigits(TClonesArray digits, AliPHOSRawFitterv0 fitter)
	157	{
	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	//--------------------------------------------------------------------------------------
	166	void AliPHOSRawDigiProducer::MakeDigits(TClonesArray digits, TClonesArray tmpDigLG, AliPHOSRawFitterv0* fitter)
	167	{
	168	//Makes the job.
	169	//TClonesArray digits, tmpDigLG and raw data fitter should be provided by calling function.
	170
	171	digits->Clear();
	172
	173	Int_t iDigit=0 ;
	174	Int_t relId[4], absId=-1, caloFlag=-1;
	175
	176	const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch.
	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() ;
	185
	186	// Clear a temporary array for LowGain digits
	187	tmpDigLG->Clear();
	188	Int_t ilgDigit=0 ;
	189
	190	//Let fitter subtract pedestals in case of ZS
	191	fitter->SetCalibData(fgCalibData) ;
	192
	193	while (fRawStream->NextDDL()) {
	194	while (fRawStream->NextChannel()) {
	195	relId[0] = 5 - fRawStream->GetModule() ; // counts from 1 to 5
	196	relId[1] = 0;
	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
	200
	201	if(caloFlag!=0 && caloFlag!=1) continue; //TRU data!
	202
	203	fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag);
	204
	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
	217	fGeom->RelToAbsNumbering(relId, absId);
	218
	219	fitter->SetNBunches(0);
	220	Int_t sigStart =0 ;
	221	Int_t sigLength=0 ;
	222	while (fRawStream->NextBunch()) { //Take the first in time bunch
	223	const UShort_t *sig = fRawStream->GetSignals();
	224	sigStart = fRawStream->GetStartTimeBin();
	225	sigLength = fRawStream->GetBunchLength();
	226	fitter->Eval(sig,sigStart,sigLength);
	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	}
	239	} // End of NextBunch()
	240
	241
	242	Double_t energy = fitter->GetEnergy() ;
	243	Double_t time = fitter->GetTime() ;
	244	if(energy<=baseLine) //in ADC channels
	245	continue ;
	246
	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
	252
	253	if(fitter->GetSignalQuality() > fSampleQualityCut && !(fitter->IsOverflow()))
	254	continue ;
	255
	256	energy = CalibrateE(energy,relId,!caloFlag) ;
	257
	258	//convert time from sample bin units to s
	259	time*=fSampleToSec ;
	260	//CalibrateT moved to Clusterizer
	261	// time = CalibrateT(time,relId,!caloFlag) ;
	262	// subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0:
	263	time -= fRawStream->GetL1Phase();
	264
	265
	266
	267	if(energy <= 0.)
	268	continue;
	269
	270	if (caloFlag == AliCaloRawStreamV3::kLowGain) {
	271	new((*tmpDigLG)[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
	272	if (sigLength>0 && fADCValuesLG!=0)
	273	static_cast<AliPHOSDigit*>(tmpDigLG->At(ilgDigit))->SetALTROSamplesLG(sigLength,fADCValuesLG);
	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);
	282	if (sigLength>0 && fADCValuesHG!=0)
	283	static_cast<AliPHOSDigit*>(digits->At(iDigit))->SetALTROSamplesHG(sigLength,fADCValuesHG);
	284	iDigit++;
	285	}
	286	} // End of NextChannel()
	287
	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() ;
	292	tmpDigLG->Sort() ;
	293	Int_t iLG = 0;
	294	Int_t nLG1 = tmpDigLG->GetEntriesFast()-1 ;
	295
	296	for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) {
	297	AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
	298	if (!digHG) continue;
	299	AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpDigLG->At(iLG)) ;
	300	while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
	301	iLG++ ;
	302	digLG = dynamic_cast<AliPHOSDigit*>(tmpDigLG->At(iLG)) ;
	303	}
	304	absId=digHG->GetId() ;
	305	fGeom->AbsToRelNumbering(absId,relId) ;
	306
	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	}
	312	}
	313	else{ //no pair - remove
	314	if(digHG->GetEnergy()<0.) //no pair, in saturation
	315	digits->RemoveAt(iDig) ;
	316	}
	317	}
	318	} // End of NextDDL()
	319
	320	CleanDigits(digits) ;
	321
	322	}
	323	//____________________________________________________________________________
	324	Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
	325	{
	326	// Convert EMC LG amplitude to HG (multipli by ~16)
	327	// Calibration parameters are taken from calibration data base
	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);
	335	}
	336	return amp ;
	337	}
	338	}
	339	return 0;
	340	}
	341	//____________________________________________________________________________
	342	Double_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	//____________________________________________________________________________
	358	void 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) ;
	382	if(fgCalibData->IsBadChannelEmc(relid[0],relid[3],relid[2])){
	383	digits->RemoveAt(i) ;
	384	}
	385	}
	386	}
	387
	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	//____________________________________________________________________________
	397	Bool_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	//____________________________________________________________________________
	405	Bool_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	//____________________________________________________________________________
	411	void 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");
	427	}