[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:
      //Very strange behaviour of electronics, but cross-checkes several times...
      if( fRawStream->GetL1Phase()<55.*1.e-9 ) //for phase=0,25,50
        time -= fRawStream->GetL1Phase();
      else //for phase 75
        time += 25.*1.e-9 ;
      
      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()) ;
	  digHG->SetLG(kTRUE) ;
	}
      }
      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
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.),
6f47f50d	55	fSampleToSec(0.),
7e88424f	56	fEmcCrystals(0),
7e88424f	57	fGeom(0),
379c5c09	58	fPulseGenerator(0),
379c5c09	59	fRawReader(0),
54ac223e	60	fRawStream(0),
	61	fADCValuesLG(0),
	62	fADCValuesHG(0)
7e88424f	63	{
7e88424f	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	//--------------------------------------------------------------------------------------
	75	AliPHOSRawDigiProducer::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	103	AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):
	104	TObject(),
	105	fEmcMinE(0.),
	106	fCpvMinE(0.),
	107	fSampleQualityCut(1.),
6f47f50d	108	fSampleToSec(0.),
7e88424f	109	fEmcCrystals(0),
7e88424f	110	fGeom(0),
379c5c09	111	fPulseGenerator(0),
379c5c09	112	fRawReader(0),
54ac223e	113	fRawStream(0),
	114	fADCValuesLG(0),
	115	fADCValuesHG(0)
	116
7e88424f	117	{
7e88424f	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	129	AliPHOSRawDigiProducer& 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 ;
c0394bfb	140	fGeom = dp.fGeom ;
39569d36	141	if(fPulseGenerator) delete fPulseGenerator ;
39569d36	142	fPulseGenerator = new AliPHOSPulseGenerator();
c0394bfb	143	return *this;
c0394bfb	144	}
379c5c09	145	//--------------------------------------------------------------------------------------
7e88424f	146	AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer()
7e88424f	147	{
379c5c09	148	// Destructor
7e88424f	149	if(fPulseGenerator) delete fPulseGenerator ;
7e88424f	150	fPulseGenerator=0 ;
379c5c09	151	delete fRawStream;
54ac223e	152	delete [] fADCValuesLG;
54ac223e	153	delete [] fADCValuesHG;
39569d36	154	}
7bc140a5	155	//--------------------------------------------------------------------------------------
379c5c09	156	void AliPHOSRawDigiProducer::MakeDigits(TClonesArray digits, AliPHOSRawFitterv0 fitter)
771123c7	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)
88fb5e50	167	{
7bc140a5	168	//Makes the job.
771123c7	169	//TClonesArray digits, tmpDigLG and raw data fitter should be provided by calling function.
7bc140a5	170
88fb5e50	171	digits->Clear();
88fb5e50	172
379c5c09	173	Int_t iDigit=0 ;
379c5c09	174	Int_t relId[4], absId=-1, caloFlag=-1;
fdd91c5a	175
c0394bfb	176	const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch.
6f47f50d	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() ;
fdd91c5a	185
771123c7	186	// Clear a temporary array for LowGain digits
771123c7	187	tmpDigLG->Clear();
c0394bfb	188	Int_t ilgDigit=0 ;
bafc1087	189
379c5c09	190	//Let fitter subtract pedestals in case of ZS
379c5c09	191	fitter->SetCalibData(fgCalibData) ;
c0394bfb	192
379c5c09	193	while (fRawStream->NextDDL()) {
379c5c09	194	while (fRawStream->NextChannel()) {
fc44a661	195	relId[0] = 5 - fRawStream->GetModule() ; // counts from 1 to 5
379c5c09	196	relId[1] = 0;
1281504a	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
f4ca5622	200
	201	if(caloFlag!=0 && caloFlag!=1) continue; //TRU data!
	202
e10f780e	203	fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag);
fc44a661	204
c4dad9d2	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
379c5c09	217	fGeom->RelToAbsNumbering(relId, absId);
f4ca5622	218
c01c5641	219	fitter->SetNBunches(0);
c2244d6d	220	Int_t sigStart =0 ;
	221	Int_t sigLength=0 ;
	222	while (fRawStream->NextBunch()) { //Take the first in time bunch
379c5c09	223	const UShort_t *sig = fRawStream->GetSignals();
54ac223e	224	sigStart = fRawStream->GetStartTimeBin();
54ac223e	225	sigLength = fRawStream->GetBunchLength();
92236b27	226	fitter->Eval(sig,sigStart,sigLength);
54ac223e	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	}
379c5c09	239	} // End of NextBunch()
70d93620	240
379c5c09	241
	242	Double_t energy = fitter->GetEnergy() ;
	243	Double_t time = fitter->GetTime() ;
	244	if(energy<=baseLine) //in ADC channels
	245	continue ;
88fb5e50	246
379c5c09	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
c0394bfb	252
379c5c09	253	if(fitter->GetSignalQuality() > fSampleQualityCut && !(fitter->IsOverflow()))
	254	continue ;
	255
817744df	256	energy = CalibrateE(energy,relId,!caloFlag) ;
6f47f50d	257
	258	//convert time from sample bin units to s
	259	time*=fSampleToSec ;
b879fc0b	260	//CalibrateT moved to Clusterizer
	261	// time = CalibrateT(time,relId,!caloFlag) ;
	262	// subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0:
4c849151	263	//Very strange behaviour of electronics, but cross-checkes several times...
	264	if( fRawStream->GetL1Phase()<55.*1.e-9 ) //for phase=0,25,50
	265	time -= fRawStream->GetL1Phase();
	266	else //for phase 75
	267	time += 25.*1.e-9 ;
379c5c09	268
	269	if(energy <= 0.)
	270	continue;
	271
	272	if (caloFlag == AliCaloRawStreamV3::kLowGain) {
771123c7	273	new((*tmpDigLG)[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
c2244d6d	274	if (sigLength>0 && fADCValuesLG!=0)
09828f52	275	static_cast<AliPHOSDigit*>(tmpDigLG->At(ilgDigit))->SetALTROSamplesLG(sigLength,fADCValuesLG);
379c5c09	276	ilgDigit++ ;
	277	}
	278	else if (caloFlag == AliCaloRawStreamV3::kHighGain) {
	279	if(fitter->IsOverflow()) //Keep this digit to replace it by Low Gain later.
	280	//If there is no LogGain it wil be removed by cut on Min E
	281	new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
	282	else
	283	new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
c2244d6d	284	if (sigLength>0 && fADCValuesHG!=0)
09828f52	285	static_cast<AliPHOSDigit*>(digits->At(iDigit))->SetALTROSamplesHG(sigLength,fADCValuesHG);
379c5c09	286	iDigit++;
	287	}
	288	} // End of NextChannel()
c0394bfb	289
379c5c09	290	//Now scan created LG and HG digits and keep only those which appeared in both lists
	291	//replace energy of HighGain digits only if there is overflow
	292	//negative energy (overflow)
	293	digits->Sort() ;
771123c7	294	tmpDigLG->Sort() ;
379c5c09	295	Int_t iLG = 0;
771123c7	296	Int_t nLG1 = tmpDigLG->GetEntriesFast()-1 ;
379c5c09	297
	298	for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) {
	299	AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
	300	if (!digHG) continue;
771123c7	301	AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpDigLG->At(iLG)) ;
379c5c09	302	while(digLG && iLG<nLG1 && digHG->GetId()> digLG->GetId()){
379c5c09	303	iLG++ ;
771123c7	304	digLG = dynamic_cast<AliPHOSDigit*>(tmpDigLG->At(iLG)) ;
379c5c09	305	}
	306	absId=digHG->GetId() ;
	307	fGeom->AbsToRelNumbering(absId,relId) ;
a7c985dc	308
379c5c09	309	if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
	310	if(digHG->GetEnergy()<0.){ //This is overflow in HG
	311	digHG->SetTime(digLG->GetTime()) ;
	312	digHG->SetEnergy(digLG->GetEnergy()) ;
f1731579	313	digHG->SetLG(kTRUE) ;
379c5c09	314	}
379c5c09	315	}
379c5c09	316	else{ //no pair - remove
fdd91c5a	317	if(digHG->GetEnergy()<0.) //no pair, in saturation
fdd91c5a	318	digits->RemoveAt(iDig) ;
88fb5e50	319	}
88fb5e50	320	}
379c5c09	321	} // End of NextDDL()
c0394bfb	322
c0394bfb	323	CleanDigits(digits) ;
379c5c09	324
c0394bfb	325	}
	326	//____________________________________________________________________________
	327	Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
	328	{
634d2178	329	// Convert EMC LG amplitude to HG (multipli by ~16)
634d2178	330	// Calibration parameters are taken from calibration data base
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	if(isLowGain){
	337	amp*= fgCalibData->GetHighLowRatioEmc(module,column,row);
88fb5e50	338	}
c0394bfb	339	return amp ;
	340	}
	341	}
	342	return 0;
	343	}
	344	//____________________________________________________________________________
	345	Double_t AliPHOSRawDigiProducer::CalibrateT(Double_t time, Int_t * relId, Bool_t /* isLowGain */)
	346	{
	347	//Calibrate time
	348	if(fgCalibData){
	349	Int_t module = relId[0];
	350	Int_t column = relId[3];
	351	Int_t row = relId[2];
	352	if(relId[1]==0) { // this is EMC
	353	time += fgCalibData->GetTimeShiftEmc(module,column,row);
	354	return time ;
	355	}
	356	}
	357
	358	return -999.;
	359	}
	360	//____________________________________________________________________________
	361	void AliPHOSRawDigiProducer::CleanDigits(TClonesArray * digits)
	362	{
	363	// remove digits with amplitudes below threshold.
	364	// remove digits in bad channels
	365	// sort digits with icreasing AbsId
	366
	367	//remove digits in bad map and below threshold
	368	Bool_t isBadMap = 0 ;
	369	if(fgCalibData->GetNumOfEmcBadChannels()){
	370	isBadMap=1 ;
	371	}
	372
	373	for(Int_t i=0; i<digits->GetEntriesFast(); i++){
	374	AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
	375	if(!digit)
	376	continue ;
	377	if ( (IsInEMC(digit) && digit->GetEnergy() < fEmcMinE) \|\|
	378	(IsInCPV(digit) && digit->GetEnergy() < fCpvMinE) ){
	379	digits->RemoveAt(i) ;
	380	continue ;
	381	}
	382	if(isBadMap){ //check bad map now
	383	Int_t relid[4] ;
	384	fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
a65517aa	385	if(fgCalibData->IsBadChannelEmc(relid[0],relid[3],relid[2])){
c0394bfb	386	digits->RemoveAt(i) ;
88fb5e50	387	}
88fb5e50	388	}
88fb5e50	389	}
88fb5e50	390
c0394bfb	391	//Compress, sort and set indexes
	392	digits->Compress() ;
	393	// digits->Sort(); already sorted earlier
	394	for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
	395	AliPHOSDigit digit = static_cast<AliPHOSDigit>( digits->At(i) ) ;
	396	digit->SetIndexInList(i) ;
	397	}
	398	}
	399	//____________________________________________________________________________
	400	Bool_t AliPHOSRawDigiProducer::IsInEMC(AliPHOSDigit * digit) const
	401	{
	402	// Tells if (true) or not (false) the digit is in a PHOS-EMC module
	403	return digit->GetId() <= fEmcCrystals ;
	404
	405	}
	406
	407	//____________________________________________________________________________
	408	Bool_t AliPHOSRawDigiProducer::IsInCPV(AliPHOSDigit * digit) const
	409	{
	410	// Tells if (true) or not (false) the digit is in a PHOS-CPV module
	411	return digit->GetId() > fEmcCrystals ;
	412	}
	413	//____________________________________________________________________________
	414	void AliPHOSRawDigiProducer::GetCalibrationParameters()
	415	{
	416	// Set calibration parameters:
	417	// if calibration database exists, they are read from database,
	418	// otherwise, reconstruction stops in the constructor of AliPHOSCalibData
	419	//
	420	// It is a user responsilibity to open CDB before reconstruction, for example:
	421	// AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
	422
	423	if (!fgCalibData){
	424	fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
	425	}
	426	if (fgCalibData->GetCalibDataEmc() == 0)
	427	AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
	428	if (fgCalibData->GetCalibDataCpv() == 0)
	429	AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
88fb5e50	430	}