X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALRawUtils.cxx;h=02a492d3ada4d7fefd42fe262c1e0651beca4d2b;hb=3b8fd9fe1c7b2d56d712b2f1086ddea59a4c4d27;hp=684b26eedc2f639da6fb9442b26ae2e950df5d59;hpb=03a3ed3ff5d505412317a6672d2890f1fbc54427;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALRawUtils.cxx b/EMCAL/AliEMCALRawUtils.cxx index 684b26eedc2..02a492d3ada 100644 --- a/EMCAL/AliEMCALRawUtils.cxx +++ b/EMCAL/AliEMCALRawUtils.cxx @@ -48,11 +48,21 @@ class AliCaloAltroMapping; class AliEMCALDigitizer; #include "AliEMCALDigit.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" + ClassImp(AliEMCALRawUtils) // Signal shape parameters -Int_t AliEMCALRawUtils::fgTimeBins = 256; // number of time bins for EMCAL +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 = 1.5E-6 ; // 15 time bins ~ 1.5 musec @@ -62,18 +72,22 @@ Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule Int_t AliEMCALRawUtils::fgPedestalValue = 32; // pedestal value for digits2raw Double_t AliEMCALRawUtils::fgFEENoise = 3.; // 3 ADC channels of noise (sampled) -AliEMCALRawUtils::AliEMCALRawUtils() +AliEMCALRawUtils::AliEMCALRawUtils(fitAlgorithm fitAlgo) : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), - fNPedSamples(0), fGeom(0), fOption("") + fNPedSamples(0), fGeom(0), fOption(""), + fRemoveBadChannels(kTRUE),fFittingAlgorithm(0),fRawAnalyzer(0) { //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 = 4; - fNPedSamples = 5; + fNoiseThreshold = 3; // 3 ADC counts is approx. noise level + fNPedSamples = 4; // less than this value => likely pedestal samples + fRemoveBadChannels = kTRUE; //Remove bad channels before fitting + SetFittingAlgorithm(fitAlgo); //Get Mapping RCU files from the AliEMCALRecParam const TObjArray* maps = AliEMCALRecParam::GetMappings(); @@ -99,9 +113,10 @@ AliEMCALRawUtils::AliEMCALRawUtils() } //____________________________________________________________________________ -AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry) +AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry, fitAlgorithm fitAlgo) : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), - fNPedSamples(0), fGeom(pGeometry), fOption("") + fNPedSamples(0), fGeom(pGeometry), fOption(""), + fRemoveBadChannels(kTRUE),fFittingAlgorithm(0),fRawAnalyzer() { // // Initialize with the given geometry - constructor required by HLT @@ -112,12 +127,16 @@ AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry) //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; - fNPedSamples = 5; - + fNoiseThreshold = 3; // 3 ADC counts is approx. noise level + fNPedSamples = 4; // less than this value => likely pedestal samples + fRemoveBadChannels = kTRUE; //Remove bad channels before fitting + SetFittingAlgorithm(fitAlgo); + + //Get Mapping RCU files from the AliEMCALRecParam const TObjArray* maps = AliEMCALRecParam::GetMappings(); if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); @@ -139,7 +158,10 @@ AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU) fNoiseThreshold(rawU.fNoiseThreshold), fNPedSamples(rawU.fNPedSamples), fGeom(rawU.fGeom), - fOption(rawU.fOption) + fOption(rawU.fOption), + fRemoveBadChannels(rawU.fRemoveBadChannels), + fFittingAlgorithm(rawU.fFittingAlgorithm), + fRawAnalyzer(rawU.fRawAnalyzer) { //copy ctor fMapping[0] = rawU.fMapping[0]; @@ -161,6 +183,9 @@ AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU) fNPedSamples = rawU.fNPedSamples; fGeom = rawU.fGeom; fOption = rawU.fOption; + fRemoveBadChannels = rawU.fRemoveBadChannels; + fFittingAlgorithm = rawU.fFittingAlgorithm; + fRawAnalyzer = rawU.fRawAnalyzer; fMapping[0] = rawU.fMapping[0]; fMapping[1] = rawU.fMapping[1]; fMapping[2] = rawU.fMapping[2]; @@ -284,7 +309,7 @@ void AliEMCALRawUtils::Digits2Raw() } //____________________________________________________________________________ -void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr) +void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap) { // convert raw data of the current event to digits @@ -301,106 +326,112 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr) AliCaloRawStreamV3 in(reader,"EMCAL",fMapping); // Select EMCAL DDL's; - reader->Select("EMCAL"); - - //Updated fitting routine from 2007 beam test takes into account - //possibility of two peaks in data and selects first one for fitting - //Also sets some of the starting parameters based on the shape of the - //given raw signal being fit - - TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); - signalF->SetParameters(10.,0.,fTau,fOrder,5.); //set all defaults once, just to be safe - signalF->SetParNames("amp","t0","tau","N","ped"); - signalF->SetParameter(2,fTau); // tau in units of time bin - signalF->SetParLimits(2,2,-1); - signalF->SetParameter(3,fOrder); // order - signalF->SetParLimits(3,2,-1); - - Int_t id = -1; - Float_t time = 0. ; - Float_t amp = 0. ; - Int_t i = 0; - Int_t startBin = 0; + reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL - //Graph to hold data we will fit (should be converted to an array - //later to speed up processing - TGraph * gSig = new TGraph(GetRawFormatTimeBins()); + // fRawAnalyzer setup + 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()) { while (in.NextChannel()) { - - //Check if the signal is high or low gain and then do the fit, - //if it is from TRU do not fit - caloFlag = in.GetCaloFlag(); - if (caloFlag != 0 && caloFlag != 1) continue; - - // There can be zero-suppression in the raw data, - // so set up the TGraph in advance - for (i=0; i < GetRawFormatTimeBins(); i++) { - gSig->SetPoint(i, i , 0); - } - - Int_t maxTime = 0; - Int_t min = 0x3ff; // init to 10-bit max - Int_t max = 0; // init to 10-bit min - int nsamples = 0; - while (in.NextBunch()) { - const UShort_t *sig = in.GetSignals(); - startBin = in.GetStartTimeBin(); - - if (((UInt_t) maxTime) < in.GetStartTimeBin()) { - maxTime = in.GetStartTimeBin(); // timebins come in reverse order - } - if (maxTime < 0 || maxTime >= GetRawFormatTimeBins()) { - AliWarning(Form("Invalid time bin %d",maxTime)); - maxTime = GetRawFormatTimeBins(); - } - nsamples += in.GetBunchLength(); - for (i = 0; i < in.GetBunchLength(); i++) { - time = startBin--; - gSig->SetPoint(time, time, (Double_t) sig[i]) ; - if (max < sig[i]) max= sig[i]; - if (min > sig[i]) min = sig[i]; - } + //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; + + //Do not fit bad channels + if(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 (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout - id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; - lowGain = in.IsLowGain(); + Float_t time = 0; + Float_t amp = 0; - gSig->Set(maxTime+1); - if ( (max - min) > fNoiseThreshold) FitRaw(gSig, signalF, amp, time) ; - - //if (caloFlag == 0 || caloFlag == 1) { // low gain or high gain - if (amp > 0 && amp < 2000) { //check both high and low end of - //result, 2000 is somewhat arbitrary - not nice with magic numbers in the code.. - AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); - - AddDigit(digitsArr, id, lowGain, (Int_t)amp, time); - } + if ( fFittingAlgorithm == kFastFit || fFittingAlgorithm == kNeuralNet || fFittingAlgorithm == kLMS || fFittingAlgorithm == kPeakFinder || fFittingAlgorithm == 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.GetTof(); + } + 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 ampEstimate = 0; + short maxADC = 0; + short timeEstimate = 0; + Float_t pedEstimate = 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; + amp = ampEstimate; + + if ( nsamples > 1 ) { // possibly something to fit + FitRaw(first, last, amp, time); + } + + if ( amp>0 && time>0 ) { // brief sanity check of fit results + + // check fit results: should be consistent with initial estimates + // more magic numbers, but very loose cuts, for now.. + // We have checked that amp and ampEstimate values are positive so division for assymmetry + // calculation should be OK/safe + Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate); + if ( (TMath::Abs(ampAsymm) > 0.1) ) { + AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations ped %f max-ped %f time %d", + amp, time, pedEstimate, ampEstimate, timeEstimate)); + + // what should do we do then? skip this channel or assign the simple estimate? + // for now just overwrite the fit results with the simple estimate + amp = ampEstimate; + time = timeEstimate; + } // asymm check + } // amp & time check + } // ampEstimate check + } // method selection + + if (amp > fNoiseThreshold) { // something to be stored + Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; + lowGain = in.IsLowGain(); - // Reset graph - for (Int_t index = 0; index < gSig->GetN(); index++) { - gSig->SetPoint(index, index, 0) ; - } - // Reset starting parameters for fit function - signalF->SetParameters(10.,0.,fTau,fOrder,5.); //reset all defaults just to be safe + // go from time-bin units to physical time fgtimetrigger + time = time * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger? - } // nsamples>0 check, some data found for this channel; not only trailer/header + 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); + // round off amplitude value to nearest integer + AddDigit(digitsArr, id, lowGain, TMath::Nint(amp), time); + } + } // end while over channel } //end while over DDL's, of input stream - - delete signalF ; - delete gSig; - + return ; } @@ -414,7 +445,6 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain // Called by Raw2Digits AliEMCALDigit *digit = 0, *tmpdigit = 0; - TIter nextdigit(digitsArr); while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { if (tmpdigit->GetId() == id) @@ -422,7 +452,7 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain } if (!digit) { // no digit existed for this tower; create one - if (lowGain) + if (lowGain && amp > fgkOverflowCut) amp = Int_t(fHighLowGainFactor * amp); Int_t idigit = digitsArr->GetEntries(); new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ; @@ -443,106 +473,157 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain } //____________________________________________________________________________ -void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_t & time) const -{ - // Fits the raw signal time distribution; from AliEMCALGetter - printf("*********************** FIT Signal\n"); - amp = time = 0. ; - Double_t ped = 0; - Int_t nPed = 0; - - for (Int_t index = 0; index < fNPedSamples; index++) { - Double_t ttime, signal; - gSig->GetPoint(index, ttime, signal) ; - if (signal > 0) { - ped += signal; - nPed++; - } - } - - if (nPed > 0) - ped /= nPed; - else { - AliWarning("Could not determine pedestal"); - ped = 10; // put some small value as first guess - } - - - Int_t maxFound = 0; - Int_t iMax = 0; - Float_t max = -1; - Float_t maxFit = gSig->GetN(); - Float_t minAfterSig = 9999; - Int_t tminAfterSig = gSig->GetN(); - Int_t nPedAfterSig = 0; - Int_t plateauWidth = 0; - Int_t plateauStart = 9999; - Float_t cut = 0.3; - - for (Int_t i=fNPedSamples; i < gSig->GetN(); i++) { - Double_t ttime, signal; - gSig->GetPoint(i, ttime, signal) ; - if (!maxFound && signal > max) { - iMax = i; - max = signal; - } - else if ( max > ped + fNoiseThreshold ) { - maxFound = 1; - minAfterSig = signal; - tminAfterSig = i; - } - if (maxFound) { - if ( signal < minAfterSig) { - minAfterSig = signal; - tminAfterSig = i; - } - if (i > tminAfterSig + 5) { // Two close peaks; end fit at minimum - maxFit = tminAfterSig; - break; +void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time) const +{ // Fits the raw signal time distribution + + //-------------------------------------------------- + //Do the fit, different fitting algorithms available + //-------------------------------------------------- + int nsamples = lastTimeBin - firstTimeBin + 1; + + switch(fFittingAlgorithm) { + case 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(timebin, timebin, fRawAnalyzer->GetReversed(timebin)); } - if ( signal < cut*max){ //stop fit at 30% amplitude(avoid the pulse shape falling edge) - maxFit = i; - break; + + 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); + + gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points + + // assign fit results + amp = signalF->GetParameter(0); + time = signalF->GetParameter(1); + + delete signalF; + + // cross-check with ParabolaFit to see if the results make sense + FitParabola(gSig, amp); // amp is possibly updated + + //printf("Std : Amp %f, time %g\n",amp, time); + delete gSig; // delete TGraph + + break; + }//kStandard Fitter + //---------------------------- + case 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) ) ); } - if ( signal < ped + fNoiseThreshold) - nPedAfterSig++; - if (nPedAfterSig >= 5) { // include 5 pedestal bins after peak - maxFit = i; - break; + + 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)); } - } - //Add check on plateau - if (signal >= fgkRawSignalOverflow - fNoiseThreshold) { - if(plateauWidth == 0) plateauStart = i; - plateauWidth++; - } - } + + 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); + + delete signalFLog; + //printf("LogFit: Amp %f, time %g\n",amp, time); + delete gSigLog; + break; + } //kLogFit + //---------------------------- + + //---------------------------- + }//switch fitting algorithms + + return; +} - if(plateauWidth > 0) { - for(int j = 0; j < plateauWidth; j++) { - //Note, have to remove the same point N times because after each - //remove, the positions of all subsequent points have shifted down - gSig->RemovePoint(plateauStart); +//__________________________________________________________________ +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 ; } } - - if ( max - ped > fNoiseThreshold ) { // else its noise - AliDebug(2,Form("Fitting max %d ped %d", max, ped)); - signalF->SetRange(0,maxFit); - - if(max-ped > 50) - signalF->SetParLimits(2,1,3); - - signalF->SetParameter(4, ped) ; - signalF->SetParameter(1, iMax); - signalF->SetParameter(0, max); - - gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points - amp = signalF->GetParameter(0); - time = signalF->GetParameter(1)*GetRawFormatTimeBinWidth() - fgTimeTrigger; + // 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 + } } + + 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) { @@ -551,9 +632,9 @@ Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) // Shape of the electronics raw reponse: // It is a semi-gaussian, 2nd order Gamma function of the general form // - // t' = (t - t0 + tau) / tau - // F = A * t**N * exp( N * ( 1 - t) ) for t >= 0 - // F = 0 for t < 0 + // 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 @@ -577,8 +658,40 @@ Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) } //__________________________________________________________________ -Bool_t AliEMCALRawUtils::RawSampledResponse( -const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const +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 ; + 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 @@ -597,9 +710,13 @@ const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const 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++) { - Double_t signal = signalF.Eval(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 @@ -608,9 +725,11 @@ const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const // March 17,09 for fast fit simulations by Alexei Pavlinov. // Get from PHOS analysis. In some sense it is open questions. - Double_t noise = gRandom->Gaus(0.,fgFEENoise); - signal += noise; - + 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 ; @@ -625,3 +744,44 @@ const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const } return lowGain ; } + +//__________________________________________________________________ +void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo) +{ + //Set fitting algorithm and initialize it if this same algorithm was not set before. + + 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 == kFastFit) { + fRawAnalyzer = new AliCaloRawAnalyzerFastFit(); + } + else if (fitAlgo == kNeuralNet) { + fRawAnalyzer = new AliCaloRawAnalyzerNN(); + } + else if (fitAlgo == kLMS) { + fRawAnalyzer = new AliCaloRawAnalyzerLMS(); + } + else if (fitAlgo == kPeakFinder) { + fRawAnalyzer = new AliCaloRawAnalyzerPeakFinder(); + } + else if (fitAlgo == kCrude) { + fRawAnalyzer = new AliCaloRawAnalyzerCrude(); + } + else { + fRawAnalyzer = new AliCaloRawAnalyzer(); + } + } + +} + +