-//____________________________________________________________________________
-void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, AliCaloCalibPedestal* pedbadmap)
-{
- // convert raw data of the current event to digits
-
- digitsArr->Clear();
-
- if (!digitsArr) {
- Error("Raw2Digits", "no digits found !");
- return;
- }
- if (!reader) {
- Error("Raw2Digits", "no raw reader found !");
- return;
- }
-
- AliCaloRawStreamV3 in(reader,"EMCAL",fMapping);
- // Select EMCAL DDL's;
- reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL
-
- //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.,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
-
- Int_t id = -1;
- Float_t time = 0. ;
- Float_t amp = 0. ;
- Float_t ped = 0. ;
- Float_t ampEstimate = 0;
- Float_t timeEstimate = 0;
- Float_t pedEstimate = 0;
- Int_t i = 0;
- Int_t startBin = 0;
-
- //Graph to hold data we will fit (should be converted to an array
- //later to speed up processing
- TGraph * gSig = new TGraph(GetRawFormatTimeBins());
-
- 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;
-
- //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;
- }
-
- // 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 , -1); // init to out-of-range values
- }
-
- Int_t maxTimeBin = 0;
- Int_t min = 0x3ff; // init to 10-bit max
- Int_t max = 0; // init to 10-bit min
- while (in.NextBunch()) {
-
- const UShort_t *sig = in.GetSignals();
- startBin = in.GetStartTimeBin();
- if (maxTimeBin < startBin) {
- maxTimeBin = startBin; // timebins come in reverse order
- }
- if (maxTimeBin < 0 || maxTimeBin >= GetRawFormatTimeBins()) {
- AliWarning(Form("Invalid time bin %d",maxTimeBin));
- maxTimeBin = GetRawFormatTimeBins();
- }
-
- for (i = 0; i < in.GetBunchLength(); i++) {
- time = startBin--;
- gSig->SetPoint((Int_t)time, time, (Double_t) sig[i]) ;
- if (max < sig[i]) max = sig[i];
- if (min > sig[i]) min = sig[i];
-
- }
- } // loop over bunches
-
- gSig->Set(maxTimeBin+1); // set actual max size of TGraph
-
- //Initialize the variables, do not keep previous values.
- // not really necessary to reset all of them (only amp and time at the moment), but better safe than sorry
- amp = -1 ;
- time = -1 ;
- ped = -1;
- ampEstimate = -1 ;
- timeEstimate = -1 ;
- pedEstimate = -1;
-
- if ( (max - min) > fNoiseThreshold) {
- switch(fFittingAlgorithm)
- {
- case kStandard:
- {
- //printf("Standard fitter \n");
- FitRaw(gSig, signalF, maxTimeBin, amp, time, ped,
- ampEstimate, timeEstimate, pedEstimate);
- break;
- }
- case kFastFit:
- {
- //printf("FastFitter \n");
- Double_t eSignal = 0;
- Double_t dAmp = amp;
- Double_t dTimeEstimate = timeEstimate;
- Double_t eTimeEstimate = 0;
- Double_t eAmp = 0;
- Double_t chi2 = 0;
-
- AliCaloFastAltroFitv0::FastFit(gSig->GetX(), gSig->GetY(), gSig->GetN(),
- eSignal, fTau,
- dAmp, eAmp, dTimeEstimate, eTimeEstimate, chi2);
- amp=dAmp;
- timeEstimate = dTimeEstimate;
- //printf("FastFitter: Amp %f, time %f, eAmp %f, eTimeEstimate %f, chi2 %f\n",amp, timeEstimate,eAmp,eTimeEstimate,chi2);
-
- break;
- }
- }
- }
-
- if ( amp>0 && amp<2000 && time>0 && time<(maxTimeBin*GetRawFormatTimeBinWidth()) ) { //check both high and low end of amplitude result, and time
- //2000 is somewhat arbitrary - not nice with magic numbers in the code..
- id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
- lowGain = in.IsLowGain();
-
- // check fit results: should be consistent with initial estimates
- // more magic numbers, but very loose cuts, for now..
- // We have checked that amp an time 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) ||
- (TMath::Abs(time - timeEstimate) > 2*GetRawFormatTimeBinWidth()) ) {
- AliDebug(2,Form("Fit results ped %f amp %f time %f not consistent with expectations ped %f max-ped %f time %d",
- ped, 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;
- }
-
- 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);
- }
-
- // Reset graph
- for (Int_t index = 0; index < gSig->GetN(); index++) {
- gSig->SetPoint(index, index, -1) ;
- }
- // Reset starting parameters for fit function
- signalF->SetParameters(10.,5.,fTau,fOrder,0.); //reset all defaults just to be safe
-
- } // end while over channel
- } //end while over DDL's, of input stream
-
- delete signalF ;
- delete gSig;
-
- return ;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff