[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 AliPHOSRawDecoder. 
//
//   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 "AliPHOSRawDecoder.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSDigit.h"
#include "AliPHOSRecoParam.h"
#include "AliPHOSCalibData.h"
#include "AliPHOSPulseGenerator.h"
#include "AliLog.h"

ClassImp(AliPHOSRawDigiProducer)

AliPHOSCalibData * AliPHOSRawDigiProducer::fgCalibData  = 0 ; 

//--------------------------------------------------------------------------------------
AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():
  TObject(),
  fEmcMinE(0.),
  fCpvMinE(0.),
  fSampleQualityCut(1.),
  fGlobalAltroOffset(0),
  fGlobalAltroThreshold(0),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0)
{
  // Default constructor
}
//--------------------------------------------------------------------------------------
AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRecoParam* recoParam):
  TObject(),
  fEmcMinE(0.),
  fCpvMinE(0.),
  fSampleQualityCut(1.),
  fGlobalAltroOffset(0),
  fGlobalAltroThreshold(0),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0)
{
  // Constructor takes paramerters from recoParam

  if(!recoParam) AliFatal("Reconstruction parameters are not set!");

  fEmcMinE = recoParam->GetEMCRawDigitThreshold();
  fCpvMinE = recoParam->GetCPVMinE();
  fSampleQualityCut = recoParam->GetEMCSampleQualityCut() ;
  fGlobalAltroOffset = recoParam->GetGlobalAltroOffset() ;
  fGlobalAltroThreshold = recoParam->GetGlobalAltroThreshold() ;

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

  fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;

  fPulseGenerator = new AliPHOSPulseGenerator();

  GetCalibrationParameters() ; 
}
//--------------------------------------------------------------------------------------                       
AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):
  TObject(),
  fEmcMinE(0.),
  fCpvMinE(0.),
  fSampleQualityCut(1.),
  fGlobalAltroOffset(0),
  fGlobalAltroThreshold(0),
  fEmcCrystals(0),
  fGeom(0),
  fPulseGenerator(0)
{                                                          
  // Copy constructor

  fEmcMinE = dp.fEmcMinE ;
  fCpvMinE = dp.fCpvMinE ;
  fSampleQualityCut = dp.fSampleQualityCut;
  fGlobalAltroOffset = dp.fGlobalAltroOffset;
  fGlobalAltroThreshold = dp.fGlobalAltroThreshold;
  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 ;
  fGlobalAltroOffset = dp.fGlobalAltroOffset ;
  fGlobalAltroThreshold = dp.fGlobalAltroThreshold ;
  fEmcCrystals = dp.fEmcCrystals ;
  fGeom = dp.fGeom ;
  if(fPulseGenerator) delete fPulseGenerator ;
  fPulseGenerator = new AliPHOSPulseGenerator();
  return  *this;
} 
//--------------------------------------------------------------------------------------                                                   
AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer()
{
  // Desctructor
  if(fPulseGenerator) delete fPulseGenerator ;
  fPulseGenerator=0 ;
}
//--------------------------------------------------------------------------------------
void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* decoder) 
{
  //Makes the job.
  //TClonesArray *digits and raw data decoder should be provided by calling function.

  digits->Clear();
 
  Int_t    iDigit   = 0 ;
  Int_t relId[4], absId =0;

  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 ;


  //Read current altro offcet from RCU
  decoder->SetAmpOffset(fGlobalAltroOffset) ;

  //Read ZS threshold from RCU
  decoder->SetAmpThreshold(fGlobalAltroThreshold) ;

  //Let decoder subtract pedestals in case of ZS
  decoder->SetCalibData(fgCalibData) ;
  
  while (decoder->NextDigit()) {

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

    //remove digits with bad shape. Decoder 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(decoder->GetSampleQuality() > fSampleQualityCut && !(decoder->IsOverflow()))
       continue ;

    Bool_t lowGainFlag = decoder->IsLowGain();

    relId[0] = decoder->GetModule();
    relId[1] = 0;
    relId[2] = decoder->GetRow();
    relId[3] = decoder->GetColumn();
    fGeom->RelToAbsNumbering(relId, absId);

    Double_t time = decoder->GetTime() ;
    time = CalibrateT(time,relId,lowGainFlag) ;

    energy = CalibrateE(energy,relId,lowGainFlag) ;

    if(energy <= 0.) 
       continue;

    if(lowGainFlag){
      new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
      ilgDigit++ ; 
    }
    else{ 
      if(decoder->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++;
    }
  }

  //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) ;                                                                                                            
    }
  }

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