// Need to add concurrent high and low-gain info in the future
// No pedestal is added to the raw signal.
//*-- Author: Marco van Leeuwen (LBL)
+//*-- Major refactoring by Per Thomas Hille
-
+//#include "AliCDBManager.h"
#include "AliEMCALRawUtils.h"
-#include "TF1.h"
-#include "TGraph.h"
-#include <TRandom.h>
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliAltroBuffer.h"
#include "AliCaloConstants.h"
#include "AliCaloRawAnalyzer.h"
#include "AliCaloRawAnalyzerFactory.h"
+#include "AliEMCALRawResponse.h"
using namespace CALO;
using namespace EMCAL;
-Double_t AliEMCALRawUtils::fgTimeTrigger = 600E-9 ; // the time of the trigger as approximately seen in the data
-Int_t AliEMCALRawUtils::fgThreshold = 1;
-Int_t AliEMCALRawUtils::fgPedestalValue = 0; // pedestal value for digits2raw, default generate ZS data
-Double_t AliEMCALRawUtils::fgFEENoise = 3.; // 3 ADC channels of noise (sampled)
-
+using std::vector;
ClassImp(AliEMCALRawUtils)
fGeom(0),
fOption(""),
fRemoveBadChannels(kFALSE),
- fFittingAlgorithm(0),
+ fFittingAlgorithm(fitAlgo),
fTimeMin(-1.),
fTimeMax(1.),
fUseFALTRO(kTRUE),
fRawAnalyzer(0),
fTriggerRawDigitMaker(0x0)
-{
+{ // ctor; set up fit algo etc
SetFittingAlgorithm(fitAlgo);
- // SetFittingAlgorithm( Algo::kLMSOffline);
-
- //Get Mapping RCU files from the AliEMCALRecParam
const TObjArray* maps = AliEMCALRecParam::GetMappings();
if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
-
- for(Int_t i = 0; i < 4; i++) {
- fMapping[i] = (AliAltroMapping*)maps->At(i);
- }
-
- //To make sure we match with the geometry in a simulation file,
- //let's try to get it first. If not, take the default geometry
+ for(Int_t i = 0; i < 4; i++)
+ {
+ fMapping[i] = (AliAltroMapping*)maps->At(i);
+ }
+
AliRunLoader *rl = AliRunLoader::Instance();
- if (rl && rl->GetAliRun()) {
+ if (rl && rl->GetAliRun())
+ {
AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
- if(emcal)fGeom = emcal->GetGeometry();
- else {
+ if(emcal)
+ {
+ fGeom = emcal->GetGeometry();
+ }
+ else
+ {
+ AliDebug(1, Form("Using default geometry in raw reco"));
+ fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
+ }
+ }
+ else
+ {
AliDebug(1, Form("Using default geometry in raw reco"));
fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
}
-
- } else {
- AliDebug(1, Form("Using default geometry in raw reco"));
- fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
- }
-
- if(!fGeom) AliFatal(Form("Could not get geometry!"));
-
- fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();
-
-}
-
-
-//____________________________________________________________________________
-AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry, Algo::fitAlgorithm fitAlgo) : //fHighLowGainFactor(16.),
- //fOrder(2),
- // fTau(2.35),
- fNoiseThreshold(3),
- fNPedSamples(4),
- fGeom(pGeometry),
- fOption(""),
- fRemoveBadChannels(kFALSE),fFittingAlgorithm(0),
- fTimeMin(-1.),fTimeMax(1.),
- fUseFALTRO(kTRUE),fRawAnalyzer(0),
- fTriggerRawDigitMaker(0x0)
-{
-
-
- // Initialize with the given geometry - constructor required by HLT
- // HLT does not use/support AliRunLoader(s) instances
- // This is a minimum intervention solution
- // Comment by MPloskon@lbl.gov
- SetFittingAlgorithm(fitAlgo);
- // SetFittingAlgorithm( Algo::kLMSOffline);
- //Get Mapping RCU files from the AliEMCALRecParam
- const TObjArray* maps = AliEMCALRecParam::GetMappings();
- if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
- for(Int_t i = 0; i < 4; i++)
- {
- fMapping[i] = (AliAltroMapping*)maps->At(i);
- }
-
if(!fGeom) AliFatal(Form("Could not get geometry!"));
- fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();
+ fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();
}
-//____________________________________________________________________________
AliEMCALRawUtils::~AliEMCALRawUtils()
{
//dtor
+ delete fRawAnalyzer;
+ delete fTriggerRawDigitMaker;
}
-//____________________________________________________________________________
void AliEMCALRawUtils::Digits2Raw()
{
// convert digits of the current event to raw data
AliRunLoader *rl = AliRunLoader::Instance();
AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
-
- // get the digits
loader->LoadDigits("EMCAL");
loader->GetEvent();
TClonesArray* digits = loader->Digits() ;
return;
}
- static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
+ static const Int_t nDDL = 20*2; // 20 SM for EMCal + DCal hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
AliAltroBuffer* buffers[nDDL];
for (Int_t i=0; i < nDDL; i++)
buffers[i] = 0;
TArrayI adcValuesHigh( TIMEBINS );
// loop over digits (assume ordered digits)
- for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++)
- {
- AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
- if(!digit)
- {
- AliFatal("NULL Digit");
- }
- else
- {
- if (digit->GetAmplitude() < fgThreshold)
- {
- continue;
- }
- //get cell indices
- Int_t nSM = 0;
- Int_t nIphi = 0;
- Int_t nIeta = 0;
- Int_t iphi = 0;
- Int_t ieta = 0;
- Int_t nModule = 0;
- fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
- fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
-
- //Check which is the RCU, 0 or 1, of the cell.
- Int_t iRCU = -111;
- //RCU0
+ for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
+ AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
+ if(!digit) {
+ AliFatal("NULL Digit");
+ } else {
+ if (digit->GetAmplitude() < AliEMCALRawResponse::GetRawFormatThreshold() ) {
+ continue;
+ }
+ //get cell indices
+ Int_t nSM = 0;
+ Int_t nIphi = 0;
+ Int_t nIeta = 0;
+ Int_t iphi = 0;
+ Int_t ieta = 0;
+ Int_t nModule = 0;
+ fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
+ fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
+
+ //Check which is the RCU, 0 or 1, of the cell.
+ Int_t iRCU = -111;
if (0<=iphi&&iphi<8) iRCU=0; // first cable row
else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half;
- //second cable row
- //RCU1
else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half;
//second cable row
else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
if (iRCU<0)
Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU);
-
+
//Which DDL?
Int_t iDDL = NRCUSPERMODULE*nSM + iRCU;
if (iDDL < 0 || iDDL >= nDDL){
Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
- }
- else{
+ } else {
if (buffers[iDDL] == 0) {
// open new file and write dummy header
TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
//iRCU=1 and even (0) SM -> RCU1A.data 1
//iRCU=0 and odd (1) SM -> RCU0C.data 2
//iRCU=1 and odd (1) SM -> RCU1C.data 3
- //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl;
buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]);
buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy;
}
buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer
// calculate the time response function
} else {
- Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ;
+ Bool_t lowgain = AliEMCALRawResponse::RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(),
+ adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ;
+
if (lowgain)
- buffers[iDDL]->WriteChannel(ieta, iphi, 0, TIMEBINS, adcValuesLow.GetArray(), fgThreshold);
+ buffers[iDDL]->WriteChannel(ieta, iphi, 0, TIMEBINS, adcValuesLow.GetArray(), AliEMCALRawResponse::GetRawFormatThreshold() );
else
- buffers[iDDL]->WriteChannel(ieta,iphi, 1, TIMEBINS, adcValuesHigh.GetArray(), fgThreshold);
+ buffers[iDDL]->WriteChannel(ieta,iphi, 1, TIMEBINS, adcValuesHigh.GetArray(), AliEMCALRawResponse::GetRawFormatThreshold() );
}
}// iDDL under the limits
}//digit exists
}
-//____________________________________________________________________________
+
void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf)
{
+ // comment
AliEMCALDigit *digit = 0, *tmpdigit = 0;
TIter nextdigit(digitsArr);
}
-//____________________________________________________________________________
void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap, TClonesArray *digitsTRG, AliEMCALTriggerData* trgData)
{
-
+ //conversion of raw data to digits
if(digitsArr) digitsArr->Clear("C");
if (!digitsArr) { Error("Raw2Digits", "no digits found !");return;}
if (!reader) {Error("Raw2Digits", "no raw reader found !");return;}
AliEMCALTriggerSTURawStream inSTU(reader);
AliCaloRawStreamV3 in(reader,"EMCAL",fMapping);
- reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL
+ reader->Select("EMCAL",0,39); // 39 = AliEMCALGeoParams::fgkLastAltroDDL
fTriggerRawDigitMaker->Reset();
fTriggerRawDigitMaker->SetIO(reader, in, inSTU, digitsTRG, trgData);
fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later
Int_t lowGain = 0;
Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref.
+ Float_t bcTimePhaseCorr = 0; // for BC-based L1 phase correction
+ Int_t bcMod4 = (reader->GetBCID() % 4); // LHC uses 40 MHz, EMCal uses 10 MHz clock
+
+ //AliCDBManager* man = AliCDBManager::Instance();
+ //Int_t runNumber = man->GetRun();
+
+ Int_t runNumber = reader->GetRunNumber();
+
+ if ((runNumber >130850 ) && (bcMod4==0 || bcMod4==1))
+ bcTimePhaseCorr = -1e-7; // subtract 100 ns for certain BC values
+
while (in.NextDDL())
{
- // fprintf(fp," TP1\n");
while (in.NextChannel())
{
- // fprintf(fp," TP2\n");
caloFlag = in.GetCaloFlag();
if (caloFlag > 2) continue; // Work with ALTRO and FALTRO
if(caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow()))
AliCaloFitResults res = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2());
if(res.GetAmp() >= fNoiseThreshold )
{
- AddDigit(digitsArr, id, lowGain, res.GetAmp(), res.GetTime(), res.GetChi2(), res.GetNdf() );
+ AddDigit(digitsArr, id, lowGain, res.GetAmp(), res.GetTime()+bcTimePhaseCorr, res.GetChi2(), res.GetNdf() );
}
}//ALTRO
else if(fUseFALTRO)
void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr)
-{
+{ // rm entries with LGnoHG (unphysical), out of time window, and too bad chi2
AliEMCALDigit *digit = 0;
Int_t n = 0;
Int_t nDigits = digitsArr->GetEntriesFast();
}
-Double_t
-AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
-{
- Double_t signal = 0.;
- Double_t tau = par[2];
- Double_t n = par[3];
- Double_t ped = par[4];
- Double_t xx = ( x[0] - par[1] + tau ) / tau ;
- if (xx <= 0)
- signal = ped ;
- else
- {
- signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
- }
- return signal ;
-}
-
-
-Bool_t AliEMCALRawUtils::RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH,
- Int_t * adcL, const Int_t keyErr) const
-{
- Bool_t lowGain = kFALSE ;
- TF1 signalF("signal", RawResponseFunction, 0, TIMEBINS, 5);
- signalF.SetParameter(0, damp) ;
- signalF.SetParameter(1, (dtime + fgTimeTrigger)/ TIMEBINWITH) ;
- signalF.SetParameter(2, TAU) ;
- signalF.SetParameter(3, ORDER);
- signalF.SetParameter(4, fgPedestalValue);
-
- Double_t signal=0.0, noise=0.0;
- for (Int_t iTime = 0; iTime < TIMEBINS; iTime++) {
- signal = signalF.Eval(iTime) ;
-
- if(keyErr>0) {
- noise = gRandom->Gaus(0.,fgFEENoise);
- signal += noise;
- }
-
- adcH[iTime] = static_cast<Int_t>(signal + 0.5) ;
- if ( adcH[iTime] > MAXBINVALUE ){ // larger than 10 bits
- adcH[iTime] = MAXBINVALUE ;
- lowGain = kTRUE ;
- }
- signal /= HGLGFACTOR;
- adcL[iTime] = static_cast<Int_t>(signal + 0.5) ;
- if ( adcL[iTime] > MAXBINVALUE ) // larger than 10 bits
- adcL[iTime] = MAXBINVALUE ;
- }
- return lowGain ;
-}
-
-
-//__________________________________________________________________
-Double_t AliEMCALRawUtils::RawResponseFunctionLog(Double_t *x, Double_t *par)
-{
- Double_t signal = 0. ;
- Double_t tau = par[2];
- Double_t n = par[3];
- Double_t xx = ( x[0] - par[1] + tau ) / tau ;
-
- if (xx < 0)
- {
- signal = par[0] - n*TMath::Log(TMath::Abs(xx)) + n * (1 - xx ) ;
- }
- else
- {
- signal = par[0] + n*TMath::Log(xx) + n * (1 - xx ) ;
- }
- return signal ;
-}
-
-
-
-//__________________________________________________________________
void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo)
-{
- // fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer( Algo::kStandard );
+{ // select which fitting algo should be used
+ delete fRawAnalyzer; // delete doesn't do anything if the pointer is 0x0
fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer( fitAlgo );
fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done, for the new methods
fRawAnalyzer->SetOverflowCut ( OVERFLOWCUT );
fRawAnalyzer->SetAmpCut(fNoiseThreshold);
fRawAnalyzer->SetFitArrayCut(fNoiseThreshold);
- // return;
}