[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.),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0),
  fRawReader(0),
  fRawStream(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.),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0),
  fRawReader(rawReader),
  fRawStream(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.),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0),
  fRawReader(0),
  fRawStream(0)
{                                                          
  // Copy constructor

  fEmcMinE = dp.fEmcMinE ;
  fCpvMinE = dp.fCpvMinE ;
  fSampleQualityCut = dp.fSampleQualityCut;
  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 ;
  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;
}
//--------------------------------------------------------------------------------------
void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawFitterv0* fitter) 
{
  //Makes the job.
  //TClonesArray *digits 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. 
  const Double_t highLowDiff=2.; //Maximal difference between High and Low channels in LG adc channels 

  //Temporary array for LowGain digits
  TClonesArray tmpLG("AliPHOSDigit",10000) ;
  Int_t ilgDigit=0 ;

  //Let fitter subtract pedestals in case of ZS
  fitter->SetCalibData(fgCalibData) ;
  
  while (fRawStream->NextDDL()) {
    while (fRawStream->NextChannel()) {
      relId[0] = fRawStream->GetModule() + 1; // 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
      fGeom->RelToAbsNumbering(relId, absId);

      Int_t nBunches = 0;
      while (fRawStream->NextBunch()) {
	nBunches++;
	if (nBunches > 1) continue;
	const UShort_t *sig = fRawStream->GetSignals();
	Int_t sigStart  = fRawStream->GetStartTimeBin();
	Int_t sigLength = fRawStream->GetBunchLength();
	fitter->SetSamples(sig,sigStart,sigLength);
      } // End of NextBunch()

      fitter->SetNBunches(nBunches);
      fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag);
      fitter->Eval();
      
      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,lowGainFlag) ;
//       time   = CalibrateT(time,relId,lowGainFlag) ;
      
      if(energy <= 0.) 
	continue;
      
      if (caloFlag == AliCaloRawStreamV3::kLowGain) {
	new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
	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);
	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() ;
    tmpLG.Sort() ;
    Int_t iLG = 0;
    Int_t nLG1 = tmpLG.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*>(tmpLG.At(iLG)) ;
      while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
	iLG++ ;
	digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.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{ //Make approximate comparison of HG and LG energies
	  Double_t de = (digHG->GetEnergy()-digLG->GetEnergy()) ; 
	  if(TMath::Abs(de)>CalibrateE(double(highLowDiff),relId,1)){ //too strong difference, remove digit
	    digits->RemoveAt(iDig) ;
	  }
	}
      }
      else{ //no pair - remove
	// temporary fix for dead LG channels
	if(relId[2]%2==1 && relId[3]%16==4) 
	  continue ;
	if(digHG->GetEnergy()>CalibrateE(double(5.),relId,1)) //One can not always find LG with Amp<5 ADC ch.
	  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
  time*=fPulseGenerator->GetRawFormatTimeTrigger() ;
  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
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"
88fb5e50	34
7bc140a5	35	// --- AliRoot header files ---
88fb5e50	36	#include "AliPHOSRawDigiProducer.h"
379c5c09	37	#include "AliPHOSRawFitterv0.h"
88fb5e50	38	#include "AliPHOSGeometry.h"
88fb5e50	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
	45	ClassImp(AliPHOSRawDigiProducer)
	46
c0394bfb	47	AliPHOSCalibData * AliPHOSRawDigiProducer::fgCalibData = 0 ;
	48
	49	//--------------------------------------------------------------------------------------
7e88424f	50	AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():
	51	TObject(),
	52	fEmcMinE(0.),
	53	fCpvMinE(0.),
	54	fSampleQualityCut(1.),
	55	fEmcCrystals(0),
	56	fGeom(0),
379c5c09	57	fPulseGenerator(0),
	58	fRawReader(0),
	59	fRawStream(0)
7e88424f	60	{
7e88424f	61	// Default constructor
70d93620	62
c0394bfb	63	fGeom=AliPHOSGeometry::GetInstance() ;
	64	if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");
	65
	66	fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
39569d36	67	fPulseGenerator = new AliPHOSPulseGenerator();
c0394bfb	68	GetCalibrationParameters() ;
8be3a30a	69
c0394bfb	70	}
379c5c09	71	//--------------------------------------------------------------------------------------
	72	AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(AliRawReader *rawReader,
	73	AliAltroMapping **mapping):
	74	TObject(),
	75	fEmcMinE(0.),
	76	fCpvMinE(0.),
	77	fSampleQualityCut(1.),
	78	fEmcCrystals(0),
	79	fGeom(0),
	80	fPulseGenerator(0),
	81	fRawReader(rawReader),
	82	fRawStream(0)
	83	{
	84	// Default constructor
	85
	86	fGeom=AliPHOSGeometry::GetInstance() ;
	87	if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");
	88
	89	fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
	90	fPulseGenerator = new AliPHOSPulseGenerator();
	91	GetCalibrationParameters() ;
	92
	93	fRawStream = new AliCaloRawStreamV3(rawReader,"PHOS",mapping);
	94
	95	}
	96	//--------------------------------------------------------------------------------------
7e88424f	97	AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):
	98	TObject(),
	99	fEmcMinE(0.),
	100	fCpvMinE(0.),
	101	fSampleQualityCut(1.),
	102	fEmcCrystals(0),
	103	fGeom(0),
379c5c09	104	fPulseGenerator(0),
	105	fRawReader(0),
	106	fRawStream(0)
7e88424f	107	{
7e88424f	108	// Copy constructor
c0394bfb	109
	110	fEmcMinE = dp.fEmcMinE ;
	111	fCpvMinE = dp.fCpvMinE ;
f78c9781	112	fSampleQualityCut = dp.fSampleQualityCut;
c0394bfb	113	fEmcCrystals = dp.fEmcCrystals ;
39569d36	114	fPulseGenerator = new AliPHOSPulseGenerator();
c0394bfb	115	fGeom = dp.fGeom ;
	116	}
	117	//--------------------------------------------------------------------------------------
7e88424f	118	AliPHOSRawDigiProducer& AliPHOSRawDigiProducer::operator= (const AliPHOSRawDigiProducer &dp)
	119	{
	120	// Assign operator
	121
c0394bfb	122	if(&dp == this) return *this;
	123
	124	fEmcMinE = dp.fEmcMinE ;
	125	fCpvMinE = dp.fCpvMinE ;
70d93620	126	fSampleQualityCut = dp.fSampleQualityCut ;
c0394bfb	127	fEmcCrystals = dp.fEmcCrystals ;
c0394bfb	128	fGeom = dp.fGeom ;
39569d36	129	if(fPulseGenerator) delete fPulseGenerator ;
39569d36	130	fPulseGenerator = new AliPHOSPulseGenerator();
c0394bfb	131	return *this;
c0394bfb	132	}
379c5c09	133	//--------------------------------------------------------------------------------------
7e88424f	134	AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer()
7e88424f	135	{
379c5c09	136	// Destructor
7e88424f	137	if(fPulseGenerator) delete fPulseGenerator ;
7e88424f	138	fPulseGenerator=0 ;
379c5c09	139	delete fRawStream;
39569d36	140	}
7bc140a5	141	//--------------------------------------------------------------------------------------
379c5c09	142	void AliPHOSRawDigiProducer::MakeDigits(TClonesArray digits, AliPHOSRawFitterv0 fitter)
88fb5e50	143	{
7bc140a5	144	//Makes the job.
379c5c09	145	//TClonesArray *digits and raw data fitter should be provided by calling function.
7bc140a5	146
88fb5e50	147	digits->Clear();
88fb5e50	148
379c5c09	149	Int_t iDigit=0 ;
379c5c09	150	Int_t relId[4], absId=-1, caloFlag=-1;
88fb5e50	151
c0394bfb	152	const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch.
a7c985dc	153	const Double_t highLowDiff=2.; //Maximal difference between High and Low channels in LG adc channels
88fb5e50	154
c0394bfb	155	//Temporary array for LowGain digits
	156	TClonesArray tmpLG("AliPHOSDigit",10000) ;
	157	Int_t ilgDigit=0 ;
bafc1087	158
379c5c09	159	//Let fitter subtract pedestals in case of ZS
379c5c09	160	fitter->SetCalibData(fgCalibData) ;
c0394bfb	161
379c5c09	162	while (fRawStream->NextDDL()) {
379c5c09	163	while (fRawStream->NextChannel()) {
1281504a	164	relId[0] = fRawStream->GetModule() + 1; // counts from 1 to 5
379c5c09	165	relId[1] = 0;
1281504a	166	relId[2] = fRawStream->GetCellX() + 1; // counts from 1 to 64
	167	relId[3] = fRawStream->GetCellZ() + 1; // counts from 1 to 56
	168	caloFlag = fRawStream->GetCaloFlag(); // 0=LG, 1=HG, 2=TRU
379c5c09	169	fGeom->RelToAbsNumbering(relId, absId);
1267d56d	170
379c5c09	171	Int_t nBunches = 0;
	172	while (fRawStream->NextBunch()) {
	173	nBunches++;
	174	if (nBunches > 1) continue;
	175	const UShort_t *sig = fRawStream->GetSignals();
	176	Int_t sigStart = fRawStream->GetStartTimeBin();
	177	Int_t sigLength = fRawStream->GetBunchLength();
	178	fitter->SetSamples(sig,sigStart,sigLength);
	179	} // End of NextBunch()
70d93620	180
379c5c09	181	fitter->SetNBunches(nBunches);
	182	fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag);
	183	fitter->Eval();
	184
	185	Double_t energy = fitter->GetEnergy() ;
	186	Double_t time = fitter->GetTime() ;
	187	if(energy<=baseLine) //in ADC channels
	188	continue ;
88fb5e50	189
379c5c09	190	//remove digits with bad shape. Fitter should calculate quality so that
	191	//in default case quality [0,1], while larger values of quality mean somehow
	192	//corrupted samples, 999 means obviously corrupted sample.
	193	//It is difficult to fit samples with overflow (even setting cut on overflow values)
	194	//because too few points are left to fit. So we do not evaluate samples with overflow
c0394bfb	195
379c5c09	196	if(fitter->GetSignalQuality() > fSampleQualityCut && !(fitter->IsOverflow()))
	197	continue ;
	198
	199	// energy = CalibrateE(energy,relId,lowGainFlag) ;
	200	// time = CalibrateT(time,relId,lowGainFlag) ;
	201
	202	if(energy <= 0.)
	203	continue;
	204
	205	if (caloFlag == AliCaloRawStreamV3::kLowGain) {
	206	new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
	207	ilgDigit++ ;
	208	}
	209	else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
	210	if(fitter->IsOverflow()) //Keep this digit to replace it by Low Gain later.
	211	//If there is no LogGain it wil be removed by cut on Min E
	212	new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
	213	else
	214	new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
	215	iDigit++;
	216	}
	217	} // End of NextChannel()
c0394bfb	218
379c5c09	219	//Now scan created LG and HG digits and keep only those which appeared in both lists
	220	//replace energy of HighGain digits only if there is overflow
	221	//negative energy (overflow)
	222	digits->Sort() ;
	223	tmpLG.Sort() ;
	224	Int_t iLG = 0;
	225	Int_t nLG1 = tmpLG.GetEntriesFast()-1 ;
	226
	227	for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) {
	228	AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
	229	if (!digHG) continue;
	230	AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
	231	while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
	232	iLG++ ;
	233	digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
	234	}
	235	absId=digHG->GetId() ;
	236	fGeom->AbsToRelNumbering(absId,relId) ;
a7c985dc	237
379c5c09	238	if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
	239	if(digHG->GetEnergy()<0.){ //This is overflow in HG
	240	digHG->SetTime(digLG->GetTime()) ;
	241	digHG->SetEnergy(digLG->GetEnergy()) ;
	242	}
	243	else{ //Make approximate comparison of HG and LG energies
	244	Double_t de = (digHG->GetEnergy()-digLG->GetEnergy()) ;
	245	if(TMath::Abs(de)>CalibrateE(double(highLowDiff),relId,1)){ //too strong difference, remove digit
	246	digits->RemoveAt(iDig) ;
	247	}
	248	}
	249	}
	250	else{ //no pair - remove
	251	// temporary fix for dead LG channels
	252	if(relId[2]%2==1 && relId[3]%16==4)
	253	continue ;
	254	if(digHG->GetEnergy()>CalibrateE(double(5.),relId,1)) //One can not always find LG with Amp<5 ADC ch.
	255	digits->RemoveAt(iDig) ;
88fb5e50	256	}
88fb5e50	257	}
379c5c09	258	} // End of NextDDL()
c0394bfb	259
c0394bfb	260	CleanDigits(digits) ;
379c5c09	261
c0394bfb	262	}
	263	//____________________________________________________________________________
	264	Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
	265	{
634d2178	266	// Convert EMC LG amplitude to HG (multipli by ~16)
634d2178	267	// Calibration parameters are taken from calibration data base
c0394bfb	268	if(fgCalibData){
	269	Int_t module = relId[0];
	270	Int_t column = relId[3];
	271	Int_t row = relId[2];
	272	if(relId[1]==0) { // this is EMC
	273	if(isLowGain){
	274	amp*= fgCalibData->GetHighLowRatioEmc(module,column,row);
88fb5e50	275	}
c0394bfb	276	return amp ;
	277	}
	278	}
	279	return 0;
	280	}
	281	//____________________________________________________________________________
	282	Double_t AliPHOSRawDigiProducer::CalibrateT(Double_t time, Int_t * relId, Bool_t /* isLowGain */)
	283	{
	284	//Calibrate time
39569d36	285	time*=fPulseGenerator->GetRawFormatTimeTrigger() ;
c0394bfb	286	if(fgCalibData){
	287	Int_t module = relId[0];
	288	Int_t column = relId[3];
	289	Int_t row = relId[2];
	290	if(relId[1]==0) { // this is EMC
	291	time += fgCalibData->GetTimeShiftEmc(module,column,row);
	292	return time ;
	293	}
	294	}
	295
	296	return -999.;
	297	}
	298	//____________________________________________________________________________
	299	void AliPHOSRawDigiProducer::CleanDigits(TClonesArray * digits)
	300	{
	301	// remove digits with amplitudes below threshold.
	302	// remove digits in bad channels
	303	// sort digits with icreasing AbsId
	304
	305	//remove digits in bad map and below threshold
	306	Bool_t isBadMap = 0 ;
	307	if(fgCalibData->GetNumOfEmcBadChannels()){
	308	isBadMap=1 ;
	309	}
	310
	311	for(Int_t i=0; i<digits->GetEntriesFast(); i++){
	312	AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
	313	if(!digit)
	314	continue ;
	315	if ( (IsInEMC(digit) && digit->GetEnergy() < fEmcMinE) \|\|
	316	(IsInCPV(digit) && digit->GetEnergy() < fCpvMinE) ){
	317	digits->RemoveAt(i) ;
	318	continue ;
	319	}
	320	if(isBadMap){ //check bad map now
	321	Int_t relid[4] ;
	322	fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
a65517aa	323	if(fgCalibData->IsBadChannelEmc(relid[0],relid[3],relid[2])){
c0394bfb	324	digits->RemoveAt(i) ;
88fb5e50	325	}
88fb5e50	326	}
88fb5e50	327	}
88fb5e50	328
c0394bfb	329	//Compress, sort and set indexes
	330	digits->Compress() ;
	331	// digits->Sort(); already sorted earlier
	332	for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
	333	AliPHOSDigit digit = static_cast<AliPHOSDigit>( digits->At(i) ) ;
	334	digit->SetIndexInList(i) ;
	335	}
	336	}
	337	//____________________________________________________________________________
	338	Bool_t AliPHOSRawDigiProducer::IsInEMC(AliPHOSDigit * digit) const
	339	{
	340	// Tells if (true) or not (false) the digit is in a PHOS-EMC module
	341	return digit->GetId() <= fEmcCrystals ;
	342
	343	}
	344
	345	//____________________________________________________________________________
	346	Bool_t AliPHOSRawDigiProducer::IsInCPV(AliPHOSDigit * digit) const
	347	{
	348	// Tells if (true) or not (false) the digit is in a PHOS-CPV module
	349	return digit->GetId() > fEmcCrystals ;
	350	}
	351	//____________________________________________________________________________
	352	void AliPHOSRawDigiProducer::GetCalibrationParameters()
	353	{
	354	// Set calibration parameters:
	355	// if calibration database exists, they are read from database,
	356	// otherwise, reconstruction stops in the constructor of AliPHOSCalibData
	357	//
	358	// It is a user responsilibity to open CDB before reconstruction, for example:
	359	// AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
	360
	361	if (!fgCalibData){
	362	fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
	363	}
	364	if (fgCalibData->GetCalibDataEmc() == 0)
	365	AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
	366	if (fgCalibData->GetCalibDataCpv() == 0)
	367	AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
88fb5e50	368	}