X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALRawUtils.cxx;h=5463e0b19a05545a2229be54c87c7ae85eec0bb7;hb=6ffa51495dc03cb0acfe60d274945e7a891bde85;hp=5e52a6a792fb595fdf08b1cc8c7e2b5bb6a3836d;hpb=168c7b3cf21f0cf36791bbfd1988a02629508615;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALRawUtils.cxx b/EMCAL/AliEMCALRawUtils.cxx index 5e52a6a792f..5463e0b19a0 100644 --- a/EMCAL/AliEMCALRawUtils.cxx +++ b/EMCAL/AliEMCALRawUtils.cxx @@ -27,219 +27,98 @@ // 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 - -#include "TF1.h" -#include "TGraph.h" -#include -class TSystem; - -class AliLog; #include "AliRun.h" #include "AliRunLoader.h" -class AliCaloAltroMapping; #include "AliAltroBuffer.h" #include "AliRawReader.h" #include "AliCaloRawStreamV3.h" #include "AliDAQ.h" - #include "AliEMCALRecParam.h" #include "AliEMCALLoader.h" #include "AliEMCALGeometry.h" -class AliEMCALDigitizer; #include "AliEMCALDigit.h" #include "AliEMCALRawDigit.h" #include "AliEMCAL.h" #include "AliCaloCalibPedestal.h" -#include "AliCaloFastAltroFitv0.h" -#include "AliCaloNeuralFit.h" #include "AliCaloBunchInfo.h" #include "AliCaloFitResults.h" -#include "AliCaloRawAnalyzerFastFit.h" -#include "AliCaloRawAnalyzerNN.h" -#include "AliCaloRawAnalyzerLMS.h" -#include "AliCaloRawAnalyzerPeakFinder.h" -#include "AliCaloRawAnalyzerCrude.h" #include "AliEMCALTriggerRawDigitMaker.h" #include "AliEMCALTriggerSTURawStream.h" #include "AliEMCALTriggerData.h" +#include "AliCaloConstants.h" +#include "AliCaloRawAnalyzer.h" +#include "AliCaloRawAnalyzerFactory.h" +#include "AliEMCALRawResponse.h" -ClassImp(AliEMCALRawUtils) - -// Signal shape parameters -Int_t AliEMCALRawUtils::fgTimeBins = 256; // number of sampling bins of the raw RO signal (we typically use 15-50; theoretical max is 1k+) -Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns -Double_t AliEMCALRawUtils::fgTimeTrigger = 600E-9 ; // the time of the trigger as approximately seen in the data +using namespace CALO; +using namespace EMCAL; -// some digitization constants -Int_t AliEMCALRawUtils::fgThreshold = 1; -Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule -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) -AliEMCALRawUtils::AliEMCALRawUtils( Algo::fitAlgorithm fitAlgo) - : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), - fNPedSamples(0), fGeom(0), fOption(""), - fRemoveBadChannels(kTRUE),fFittingAlgorithm(0), - fTimeMin(-1.),fTimeMax(1.), - fUseFALTRO(kFALSE),fRawAnalyzer(0), - fTriggerRawDigitMaker(0x0) -{ - //These are default parameters. - //Can be re-set from without with setter functions - //Already set in the OCDB and passed via setter in the AliEMCALReconstructor - fHighLowGainFactor = 16. ; // Adjusted for a low gain range of 82 GeV (10 bits) - fOrder = 2; // Order of gamma fn - fTau = 2.35; // in units of timebin, from CERN 2007 testbeam - fNoiseThreshold = 3; // 3 ADC counts is approx. noise level - fNPedSamples = 4; // Less than this value => likely pedestal samples - fRemoveBadChannels = kFALSE; // Do not remove bad channels before fitting - fUseFALTRO = kTRUE; // Get the trigger FALTRO information and pass it to digits. +AliEMCALRawUtils::AliEMCALRawUtils( Algo::fitAlgorithm fitAlgo) : fNoiseThreshold(3), + fNPedSamples(4), + fGeom(0), + fOption(""), + fRemoveBadChannels(kFALSE), + fFittingAlgorithm(fitAlgo), + fTimeMin(-1.), + fTimeMax(1.), + fUseFALTRO(kTRUE), + fRawAnalyzer(0), + fTriggerRawDigitMaker(0x0) +{ // ctor; set up fit algo etc SetFittingAlgorithm(fitAlgo); - - //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(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(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), - fNPedSamples(0), fGeom(pGeometry), fOption(""), - fRemoveBadChannels(kTRUE),fFittingAlgorithm(0), - fTimeMin(-1.),fTimeMax(1.), - fUseFALTRO(kFALSE),fRawAnalyzer(), - 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 - // - //These are default parameters. - //Can be re-set from without with setter functions - //Already set in the OCDB and passed via setter in the AliEMCALReconstructor - fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) - fOrder = 2; // order of gamma fn - fTau = 2.35; // in units of timebin, from CERN 2007 testbeam - fNoiseThreshold = 3; // 3 ADC counts is approx. noise level - fNPedSamples = 4; // Less than this value => likely pedestal samples - fRemoveBadChannels = kFALSE; // Do not remove bad channels before fitting - fUseFALTRO = kTRUE; // Get the trigger FALTRO information and pass it to digits. - SetFittingAlgorithm(fitAlgo); - - //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(); -} - -//____________________________________________________________________________ -AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU) - : TObject(), - fHighLowGainFactor(rawU.fHighLowGainFactor), - fOrder(rawU.fOrder), - fTau(rawU.fTau), - fNoiseThreshold(rawU.fNoiseThreshold), - fNPedSamples(rawU.fNPedSamples), - fGeom(rawU.fGeom), - fOption(rawU.fOption), - fRemoveBadChannels(rawU.fRemoveBadChannels), - fFittingAlgorithm(rawU.fFittingAlgorithm), - fTimeMin(rawU.fTimeMin),fTimeMax(rawU.fTimeMax), - fUseFALTRO(rawU.fUseFALTRO), - fRawAnalyzer(rawU.fRawAnalyzer), - fTriggerRawDigitMaker(rawU.fTriggerRawDigitMaker) +AliEMCALRawUtils::~AliEMCALRawUtils() { - //copy ctor - fMapping[0] = rawU.fMapping[0]; - fMapping[1] = rawU.fMapping[1]; - fMapping[2] = rawU.fMapping[2]; - fMapping[3] = rawU.fMapping[3]; -} - -//____________________________________________________________________________ -AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU) -{ - //assignment operator - - if(this != &rawU) { - fHighLowGainFactor = rawU.fHighLowGainFactor; - fOrder = rawU.fOrder; - fTau = rawU.fTau; - fNoiseThreshold = rawU.fNoiseThreshold; - fNPedSamples = rawU.fNPedSamples; - fGeom = rawU.fGeom; - fOption = rawU.fOption; - fRemoveBadChannels = rawU.fRemoveBadChannels; - fFittingAlgorithm = rawU.fFittingAlgorithm; - fTimeMin = rawU.fTimeMin; - fTimeMax = rawU.fTimeMax; - fUseFALTRO = rawU.fUseFALTRO; - fRawAnalyzer = rawU.fRawAnalyzer; - fMapping[0] = rawU.fMapping[0]; - fMapping[1] = rawU.fMapping[1]; - fMapping[2] = rawU.fMapping[2]; - fMapping[3] = rawU.fMapping[3]; - fTriggerRawDigitMaker = rawU.fTriggerRawDigitMaker; - } - - return *this; - -} - -//____________________________________________________________________________ -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(rl->GetDetectorLoader("EMCAL")); - - // get the digits loader->LoadDigits("EMCAL"); loader->GetEvent(); TClonesArray* digits = loader->Digits() ; @@ -254,79 +133,82 @@ void AliEMCALRawUtils::Digits2Raw() for (Int_t i=0; i < nDDL; i++) buffers[i] = 0; - TArrayI adcValuesLow(fgTimeBins); - TArrayI adcValuesHigh(fgTimeBins); + TArrayI adcValuesLow( TIMEBINS ); + TArrayI adcValuesHigh( TIMEBINS ); // loop over digits (assume ordered digits) - for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { - AliEMCALDigit* digit = dynamic_cast(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 - 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 (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same + for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) + { + AliEMCALDigit* digit = dynamic_cast(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) ; - if (iRCU<0) - Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU); + //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; + 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 (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same + + if (iRCU<0) + Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU); - //Which DDL? - Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU; - if (iDDL < 0 || iDDL >= nDDL){ - Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); - } - else{ - if (buffers[iDDL] == 0) { - // open new file and write dummy header - TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); - //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C + //Which DDL? + Int_t iDDL = NRCUSPERMODULE*nSM + iRCU; + if (iDDL < 0 || iDDL >= nDDL){ + Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); + } + else{ + if (buffers[iDDL] == 0) + { + // open new file and write dummy header + TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); + //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C Int_t iRCUside=iRCU+(nSM%2)*2; //iRCU=0 and even (0) SM -> RCU0A.data 0 //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 "<WriteDataHeader(kTRUE, kFALSE); //Dummy; } // out of time range signal (?) - if (digit->GetTimeR() > GetRawFormatTimeMax() ) { + if (digit->GetTimeR() > TIMEBINMAX ) { AliInfo("Signal is out of time range.\n"); buffers[iDDL]->FillBuffer((Int_t)digit->GetAmplitude()); - buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin + buffers[iDDL]->FillBuffer( TIMEBINS ); // time bin buffers[iDDL]->FillBuffer(3); // bunch length 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()) ; - if (lowgain) - buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow.GetArray(), fgThreshold); + Bool_t lowgain = AliEMCALRawResponse::RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(), + adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; + + if (lowgain) + buffers[iDDL]->WriteChannel(ieta, iphi, 0, TIMEBINS, adcValuesLow.GetArray(), AliEMCALRawResponse::GetRawFormatThreshold() ); else - buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh.GetArray(), fgThreshold); + buffers[iDDL]->WriteChannel(ieta,iphi, 1, TIMEBINS, adcValuesHigh.GetArray(), AliEMCALRawResponse::GetRawFormatThreshold() ); } }// iDDL under the limits }//digit exists @@ -344,252 +226,147 @@ void AliEMCALRawUtils::Digits2Raw() loader->UnloadDigits(); } -//____________________________________________________________________________ -void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap, TClonesArray *digitsTRG, AliEMCALTriggerData* trgData) + + +void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf) { - // convert raw data of the current event to digits - - if(digitsArr) digitsArr->Clear("C"); + // comment + AliEMCALDigit *digit = 0, *tmpdigit = 0; + TIter nextdigit(digitsArr); + + while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) + { + if (tmpdigit->GetId() == id) digit = tmpdigit; + } - if (!digitsArr) { - Error("Raw2Digits", "no digits found !"); - return; - } - if (!reader) { - Error("Raw2Digits", "no raw reader found !"); - return; - } + if (!digit) { // no digit existed for this tower; create one + Int_t type = AliEMCALDigit::kHG; // use enum in AliEMCALDigit + if (lowGain) + { + amp *= HGLGFACTOR; + type = AliEMCALDigit::kLGnoHG; + } + + Int_t idigit = digitsArr->GetEntries(); + new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit, chi2, ndf); + AliDebug(2,Form("Add digit Id %d for the first time, type %d", id, type)); + }//digit added first time + else + { // a digit already exists, check range + // (use high gain if signal < cut value, otherwise low gain) + if (lowGain) + { // new digit is low gain + if (digit->GetAmplitude() > OVERFLOWCUT ) + { // use if previously stored (HG) digit is out of range + digit->SetAmplitude( HGLGFACTOR * amp); + digit->SetTime(time); + digit->SetType(AliEMCALDigit::kLG); + AliDebug(2,Form("Add LG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType())); + } + }//new low gain digit + else { // new digit is high gain + if (amp < OVERFLOWCUT ) + { // new digit is high gain; use if not out of range + digit->SetAmplitude(amp); + digit->SetTime(time); + digit->SetType(AliEMCALDigit::kHG); + AliDebug(2,Form("Add HG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType())); + } + else + { // HG out of range, just change flag value to show that HG did exist + digit->SetType(AliEMCALDigit::kLG); + AliDebug(2,Form("Change LG digit to HG, ID %d, type %d", digit->GetId(), digit->GetType())); + } + }//new high gain digit + }//digit existed replace it +} + + +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); - - // Select EMCAL DDL's; reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL - fTriggerRawDigitMaker->Reset(); fTriggerRawDigitMaker->SetIO(reader, in, inSTU, digitsTRG, trgData); - - // fRawAnalyzer setup - fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done, for the new methods - fRawAnalyzer->SetOverflowCut(fgkOverflowCut); - fRawAnalyzer->SetAmpCut(fNoiseThreshold); - fRawAnalyzer->SetFitArrayCut(fNoiseThreshold); fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later - - // channel info parameters + Int_t lowGain = 0; Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref. - // start loop over input stream - while (in.NextDDL()) { - - // if ( in.GetDDLNumber() != 0 && in.GetDDLNumber() != 2 ) continue; - - while (in.NextChannel()) { - - //Check if the signal is high or low gain and then do the fit, - //if it is from TRU or LEDMon do not fit - caloFlag = in.GetCaloFlag(); - // if (caloFlag != 0 && caloFlag != 1) continue; - if (caloFlag > 2) continue; // Work with ALTRO and FALTRO - - //Do not fit bad channels of ALTRO - if(caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) { - //printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow()); - continue; - } - - vector bunchlist; - while (in.NextBunch()) { - bunchlist.push_back( AliCaloBunchInfo(in.GetStartTimeBin(), in.GetBunchLength(), in.GetSignals() ) ); - } // loop over bunches - - if ( caloFlag < 2 ) - { // ALTRO - Float_t time = 0; - Float_t amp = 0; - short timeEstimate = 0; - Float_t ampEstimate = 0; - Bool_t fitDone = kFALSE; - Float_t chi2 = 0; - Int_t ndf = 0; - - if ( fFittingAlgorithm == Algo::kFastFit || fFittingAlgorithm == Algo::kNeuralNet || - fFittingAlgorithm == Algo::kLMS || fFittingAlgorithm == Algo::kPeakFinder || - fFittingAlgorithm == Algo::kCrude) { - // all functionality to determine amp and time etc is encapsulated inside the Evaluate call for these methods - AliCaloFitResults fitResults = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2()); - - amp = fitResults.GetAmp(); - time = fitResults.GetTime(); - timeEstimate = fitResults.GetMaxTimebin(); - ampEstimate = fitResults.GetMaxSig(); - chi2 = fitResults.GetChi2(); - ndf = fitResults.GetNdf(); - if (fitResults.GetStatus() == Ret::kFitPar) { - fitDone = kTRUE; - } - } - else { // for the other methods we for now use the functionality of - // AliCaloRawAnalyzer as well, to select samples and prepare for fits, - // if it looks like there is something to fit - - // parameters init. - Float_t pedEstimate = 0; - short maxADC = 0; - Int_t first = 0; - Int_t last = 0; - Int_t bunchIndex = 0; - // - // The PreFitEvaluateSamples + later call to FitRaw will hopefully - // be replaced by a single Evaluate call or so soon, like for the other - // methods, but this should be good enough for evaluation of - // the methods for now (Jan. 2010) - // - int nsamples = fRawAnalyzer->PreFitEvaluateSamples( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2(), bunchIndex, ampEstimate, maxADC, timeEstimate, pedEstimate, first, last); - - if (ampEstimate >= fNoiseThreshold) { // something worth looking at - - time = timeEstimate; // maxrev in AliCaloRawAnalyzer speak; comes with an offset w.r.t. real timebin - Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1); - amp = ampEstimate; - - if ( nsamples > 1 && maxADCGetBCID() % 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 ( fitDone ) { // brief sanity check of fit results - Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate); - Float_t timeDiff = time - timeEstimate; - if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) ) { - // AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations amp %f time %d", amp, time, ampEstimate, timeEstimate)); - - // for now just overwrite the fit results with the simple/initial estimate - amp = ampEstimate; - time = timeEstimate; - fitDone = kFALSE; - } - } // fitDone - - if (amp >= fNoiseThreshold) { // something to be stored - if ( ! fitDone) { // smear ADC with +- 0.5 uniform (avoid discrete effects) - amp += (0.5 - gRandom->Rndm()); // Rndm generates a number in ]0,1] - } - - Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; - lowGain = in.IsLowGain(); - - // go from time-bin units to physical time fgtimetrigger - time = time * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger? - // subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0: - time -= in.GetL1Phase(); - - AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); - // printf("Added tower: SM %d, row %d, column %d, amp %3.2f\n",in.GetModule(), in.GetRow(), in.GetColumn(),amp); + if ((runNumber >130850 ) && (bcMod4==0 || bcMod4==1)) + bcTimePhaseCorr = -1e-7; // subtract 100 ns for certain BC values - AddDigit(digitsArr, id, lowGain, amp, time, chi2, ndf); - } - - }//ALTRO - else if(fUseFALTRO) - {// Fake ALTRO - fTriggerRawDigitMaker->Add( bunchlist ); - }//Fake ALTRO - } // end while over channel - } //end while over DDL's, of input stream - + while (in.NextDDL()) + { + while (in.NextChannel()) + { + caloFlag = in.GetCaloFlag(); + if (caloFlag > 2) continue; // Work with ALTRO and FALTRO + if(caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) + { + continue; + } + vector bunchlist; + while (in.NextBunch()) + { + bunchlist.push_back( AliCaloBunchInfo(in.GetStartTimeBin(), in.GetBunchLength(), in.GetSignals() ) ); + } + if (bunchlist.size() == 0) continue; + if ( caloFlag < 2 ) + { // ALTRO + Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; + lowGain = in.IsLowGain(); + fRawAnalyzer->SetL1Phase( in.GetL1Phase() ); + AliCaloFitResults res = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2()); + if(res.GetAmp() >= fNoiseThreshold ) + { + AddDigit(digitsArr, id, lowGain, res.GetAmp(), res.GetTime()+bcTimePhaseCorr, res.GetChi2(), res.GetNdf() ); + } + }//ALTRO + else if(fUseFALTRO) + {// Fake ALTRO + fTriggerRawDigitMaker->Add( bunchlist ); + }//Fake ALTRO + } // end while over channel + } //end while over DDL's, of input stream fTriggerRawDigitMaker->PostProcess(); - TrimDigits(digitsArr); - - return ; } -//____________________________________________________________________________ -void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf) { - // - // Add a new digit. - // This routine checks whether a digit exists already for this tower - // and then decides whether to use the high or low gain info - // - // Called by Raw2Digits - - AliEMCALDigit *digit = 0, *tmpdigit = 0; - TIter nextdigit(digitsArr); - while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { - if (tmpdigit->GetId() == id) digit = tmpdigit; - } - if (!digit) { // no digit existed for this tower; create one - Int_t type = AliEMCALDigit::kHG; // use enum in AliEMCALDigit - if (lowGain) { - amp *= fHighLowGainFactor; - type = AliEMCALDigit::kLGnoHG; - } - Int_t idigit = digitsArr->GetEntries(); - new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit, chi2, ndf); - AliDebug(2,Form("Add digit Id %d for the first time, type %d", id, type)); - }//digit added first time - else { // a digit already exists, check range - // (use high gain if signal < cut value, otherwise low gain) - if (lowGain) { // new digit is low gain - if (digit->GetAmplitude() > fgkOverflowCut) { // use if previously stored (HG) digit is out of range - digit->SetAmplitude(fHighLowGainFactor * amp); - digit->SetTime(time); - digit->SetType(AliEMCALDigit::kLG); - AliDebug(2,Form("Add LG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType())); - } - }//new low gain digit - else { // new digit is high gain - if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range - digit->SetAmplitude(amp); - digit->SetTime(time); - digit->SetType(AliEMCALDigit::kHG); - AliDebug(2,Form("Add HG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType())); - } - else { // HG out of range, just change flag value to show that HG did exist - digit->SetType(AliEMCALDigit::kLG); - AliDebug(2,Form("Change LG digit to HG, ID %d, type %d", digit->GetId(), digit->GetType())); - } - }//new high gain digit - }//digit existed replace it - -} - -//____________________________________________________________________________ void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr) -{ - // Remove digits with only low gain and large time - +{ // 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(); TIter nextdigit(digitsArr); while ((digit = (AliEMCALDigit*) nextdigit())) { - - //Check if only LG existed, remove if so if (digit->GetType() == AliEMCALDigit::kLGnoHG) { AliDebug(1,Form("Remove digit with id %d, LGnoHG",digit->GetId())); digitsArr->Remove(digit); } - //Check if time is too large or too small, remove if so else if(fTimeMin > digit->GetTime() || fTimeMax < digit->GetTime()) { digitsArr->Remove(digit); AliDebug(1,Form("Remove digit with id %d, Bad Time %e",digit->GetId(), digit->GetTime())); } - // Check if Chi2 is undefined else if (0 > digit->GetChi2()) { digitsArr->Remove(digit); AliDebug(1,Form("Remove digit with id %d, Bad Chi2 %e",digit->GetId(), digit->GetChi2())); } - //Good digit, just reassign the index of the digit in case there was a previous removal else { digit->SetIndexInList(n); n++; @@ -598,366 +375,18 @@ void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr) digitsArr->Compress(); AliDebug(1,Form("N Digits before trimming : %d; after array compression %d",nDigits,digitsArr->GetEntriesFast())); - -} - -//____________________________________________________________________________ -void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const -{ // Fits the raw signal time distribution - - //-------------------------------------------------- - //Do the fit, different fitting algorithms available - //-------------------------------------------------- - int nsamples = lastTimeBin - firstTimeBin + 1; - fitDone = kFALSE; - - switch(fFittingAlgorithm) { - case Algo::kStandard: - { - if (nsamples < 3) { return; } // nothing much to fit - //printf("Standard fitter \n"); - - // Create Graph to hold data we will fit - TGraph *gSig = new TGraph( nsamples); - for (int i=0; iSetPoint(i, timebin, fRawAnalyzer->GetReversed(timebin)); - } - - TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); - signalF->SetParameters(10.,5.,fTau,fOrder,0.); //set all defaults once, just to be safe - signalF->SetParNames("amp","t0","tau","N","ped"); - signalF->FixParameter(2,fTau); // tau in units of time bin - signalF->FixParameter(3,fOrder); // order - signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here - signalF->SetParameter(1, time); - signalF->SetParameter(0, amp); - // set rather loose parameter limits - signalF->SetParLimits(0, 0.5*amp, 2*amp ); - signalF->SetParLimits(1, time - 4, time + 4); - - try { - gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points - // assign fit results - amp = signalF->GetParameter(0); - time = signalF->GetParameter(1); - chi2 = signalF->GetChisquare(); - fitDone = kTRUE; - } - catch (const std::exception & e) { - AliError( Form("TGraph Fit exception %s", e.what()) ); - // stay with default amp and time in case of exception, i.e. no special action required - fitDone = kFALSE; - } - delete signalF; - - //printf("Std : Amp %f, time %g\n",amp, time); - delete gSig; // delete TGraph - - break; - }//kStandard Fitter - //---------------------------- - case Algo::kLogFit: - { - if (nsamples < 3) { return; } // nothing much to fit - //printf("LogFit \n"); - - // Create Graph to hold data we will fit - TGraph *gSigLog = new TGraph( nsamples); - for (int i=0; iSetPoint(timebin, timebin, TMath::Log(fRawAnalyzer->GetReversed(timebin) ) ); - } - - TF1 * signalFLog = new TF1("signalLog", RawResponseFunctionLog, 0, GetRawFormatTimeBins(), 5); - signalFLog->SetParameters(2.3, 5.,fTau,fOrder,0.); //set all defaults once, just to be safe - signalFLog->SetParNames("amplog","t0","tau","N","ped"); - signalFLog->FixParameter(2,fTau); // tau in units of time bin - signalFLog->FixParameter(3,fOrder); // order - signalFLog->FixParameter(4, 0); // pedestal should be subtracted when we get here - signalFLog->SetParameter(1, time); - if (amp>=1) { - signalFLog->SetParameter(0, TMath::Log(amp)); - } - - gSigLog->Fit(signalFLog, "QROW"); // Note option 'W': equal errors on all points - - // assign fit results - Double_t amplog = signalFLog->GetParameter(0); //Not Amp, but Log of Amp - amp = TMath::Exp(amplog); - time = signalFLog->GetParameter(1); - fitDone = kTRUE; - - delete signalFLog; - //printf("LogFit: Amp %f, time %g\n",amp, time); - delete gSigLog; - break; - } //kLogFit - //---------------------------- - - //---------------------------- - }//switch fitting algorithms - - return; -} - -//__________________________________________________________________ -void AliEMCALRawUtils::FitParabola(const TGraph *gSig, Float_t & amp) const -{ - //BEG YS alternative methods to calculate the amplitude - Double_t * ymx = gSig->GetX() ; - Double_t * ymy = gSig->GetY() ; - const Int_t kN = 3 ; - Double_t ymMaxX[kN] = {0., 0., 0.} ; - Double_t ymMaxY[kN] = {0., 0., 0.} ; - Double_t ymax = 0. ; - // find the maximum amplitude - Int_t ymiMax = 0 ; - for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) { - if (ymy[ymi] > ymMaxY[0] ) { - ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude - ymMaxX[0] = ymx[ymi] ; - ymiMax = ymi ; - } - } - // find the maximum by fitting a parabola through the max and the two adjacent samples - if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) { - ymMaxY[1] = ymy[ymiMax+1] ; - ymMaxY[2] = ymy[ymiMax-1] ; - ymMaxX[1] = ymx[ymiMax+1] ; - ymMaxX[2] = ymx[ymiMax-1] ; - if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) { - //fit a parabola through the 3 points y= a+bx+x*x*x - Double_t sy = 0 ; - Double_t sx = 0 ; - Double_t sx2 = 0 ; - Double_t sx3 = 0 ; - Double_t sx4 = 0 ; - Double_t sxy = 0 ; - Double_t sx2y = 0 ; - for (Int_t i = 0; i < kN ; i++) { - sy += ymMaxY[i] ; - sx += ymMaxX[i] ; - sx2 += ymMaxX[i]*ymMaxX[i] ; - sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; - sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; - sxy += ymMaxX[i]*ymMaxY[i] ; - sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ; - } - Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2); - Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ; - Double_t c = cN / cD ; - Double_t b = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ; - Double_t a = (sy-b*sx-c*sx2)/kN ; - Double_t xmax = -b/(2*c) ; - ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude - amp = ymax; - } - } - - Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ; - if (diff > 0.1) - amp = ymMaxY[0] ; - //printf("Yves : Amp %f, time %g\n",amp, time); - //END YS - return; -} - -//__________________________________________________________________ -Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) -{ - // Matches version used in 2007 beam test - // - // Shape of the electronics raw reponse: - // It is a semi-gaussian, 2nd order Gamma function of the general form - // - // xx = (t - t0 + tau) / tau [xx is just a convenient help variable] - // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0 - // F = 0 for xx < 0 - // - // parameters: - // A: par[0] // Amplitude = peak value - // t0: par[1] - // tau: par[2] - // N: par[3] - // ped: par[4] - // - 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 ; -} - -//__________________________________________________________________ -Double_t AliEMCALRawUtils::RawResponseFunctionLog(Double_t *x, Double_t *par) -{ - // Matches version used in 2007 beam test - // - // Shape of the electronics raw reponse: - // It is a semi-gaussian, 2nd order Gamma function of the general form - // - // xx = (t - t0 + tau) / tau [xx is just a convenient help variable] - // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0 - // F = 0 for xx < 0 - // - // parameters: - // Log[A]: par[0] // Amplitude = peak value - // t0: par[1] - // tau: par[2] - // N: par[3] - // ped: par[4] - // - Double_t signal = 0. ; - Double_t tau = par[2]; - Double_t n = par[3]; - //Double_t ped = par[4]; // not used - 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 ; -} - -//__________________________________________________________________ -Bool_t AliEMCALRawUtils::RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL, const Int_t keyErr) const -{ - // for a start time dtime and an amplitude damp given by digit, - // calculates the raw sampled response AliEMCAL::RawResponseFunction - Bool_t lowGain = kFALSE ; - - // A: par[0] // Amplitude = peak value - // t0: par[1] - // tau: par[2] - // N: par[3] - // ped: par[4] - - TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); - signalF.SetParameter(0, damp) ; - signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ; - signalF.SetParameter(2, fTau) ; - signalF.SetParameter(3, fOrder); - signalF.SetParameter(4, fgPedestalValue); - - Double_t signal=0.0, noise=0.0; - for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { - signal = signalF.Eval(iTime) ; - // Next lines commeted for the moment but in principle it is not necessary to add - // extra noise since noise already added at the digits level. - - //According to Terry Awes, 13-Apr-2008 - //add gaussian noise in quadrature to each sample - //Double_t noise = gRandom->Gaus(0.,fgFEENoise); - //signal = sqrt(signal*signal + noise*noise); - - // March 17,09 for fast fit simulations by Alexei Pavlinov. - // Get from PHOS analysis. In some sense it is open questions. - if(keyErr>0) { - noise = gRandom->Gaus(0.,fgFEENoise); - signal += noise; - } - - adcH[iTime] = static_cast(signal + 0.5) ; - if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits - adcH[iTime] = fgkRawSignalOverflow ; - lowGain = kTRUE ; - } - - signal /= fHighLowGainFactor; - - adcL[iTime] = static_cast(signal + 0.5) ; - if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits - adcL[iTime] = fgkRawSignalOverflow ; - - } - - return lowGain ; } -//__________________________________________________________________ -void AliEMCALRawUtils::CalculateChi2(const Double_t* t, const Double_t* y, const Int_t nPoints, -const Double_t sig, const Double_t tau, const Double_t amp, const Double_t t0, Double_t &chi2) -{ - // Input: - // t[] - array of time bins - // y[] - array of amplitudes after pedestal subtractions; - // nPoints - number of points - // sig - error of amplitude measurement (one value for all channels) - // if sig<0 that mean sig=1. - // tau - filter time response (in timebin units) - // amp - amplitude at t0; - // t0 - time of max amplitude; - // Output: - // chi2 - chi2 - // ndf = nPoints - 2 when tau fixed - // ndf = nPoints - 3 when tau free - static Double_t par[5]={0.0, 0.0, 0.0, 2.0, 0.0}; - - par[0] = amp; - par[1] = t0; - par[2] = tau; - // par[3]=n=2.; par[4]=ped=0.0 - Double_t dy = 0.0, x = 0.0, f=0.0; - for(Int_t i=0; i f %f : dy %f \n", i, y[i], f, dy); - } - if(sig>0.0) chi2 /= (sig*sig); -} - -//__________________________________________________________________ void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo) -{ - //Set fitting algorithm and initialize it if this same algorithm was not set before. - //printf("**** Set Algorithm , number %d ****\n",fitAlgo); - - if(fitAlgo == fFittingAlgorithm && fRawAnalyzer) { - //Do nothing, this same algorithm already set before. - //printf("**** Algorithm already set before, number %d, %s ****\n",fitAlgo, fRawAnalyzer->GetName()); - return; - } - //Initialize the requested algorithm - if(fitAlgo != fFittingAlgorithm || !fRawAnalyzer) { - //printf("**** Init Algorithm , number %d ****\n",fitAlgo); - - fFittingAlgorithm = fitAlgo; - if (fRawAnalyzer) delete fRawAnalyzer; // delete prev. analyzer if existed. - - if (fitAlgo == Algo::kFastFit) { - fRawAnalyzer = new AliCaloRawAnalyzerFastFit(); - } - else if (fitAlgo == Algo::kNeuralNet) { - fRawAnalyzer = new AliCaloRawAnalyzerNN(); - } - else if (fitAlgo == Algo::kLMS) { - fRawAnalyzer = new AliCaloRawAnalyzerLMS(); - } - else if (fitAlgo == Algo::kPeakFinder) { - fRawAnalyzer = new AliCaloRawAnalyzerPeakFinder(); - } - else if (fitAlgo == Algo::kCrude) { - fRawAnalyzer = new AliCaloRawAnalyzerCrude(); - } - else { - // fRawAnalyzer = new AliCaloRawAnalyzer(); - fRawAnalyzer = 0; - } - } - +{ // 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); } +