- //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(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) {
- FitRaw(gSig, signalF, maxTimeBin, amp, time, ped,
- ampEstimate, timeEstimate, pedEstimate);
- }
-
- 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 ;
-}
-
-//____________________________________________________________________________
-void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) {
- //
- // 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
-