EMCAL/AliEMCALRawUtils.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //_________________________________________________________________________
  19 //  Utility Class for handling Raw data
  20 //  Does all transitions from Digits to Raw and vice versa,
  21 //  for simu and reconstruction
  22 //
  23 //  Note: the current version is still simplified. Only
  24 //    one raw signal per digit is generated; either high-gain or low-gain
  25 //    Need to add concurrent high and low-gain info in the future
  26 //    No pedestal is added to the raw signal.
  27 //*-- Author: Marco van Leeuwen (LBL)
  28
  29 #include "AliEMCALRawUtils.h"
  30
  31 #include "TF1.h"
  32 #include "TGraph.h"
  33 class TSystem;
  34
  35 class AliLog;
  36 #include "AliRun.h"
  37 #include "AliRunLoader.h"
  38 class AliCaloAltroMapping;
  39 #include "AliAltroBuffer.h"
  40 #include "AliRawReader.h"
  41 #include "AliCaloRawStreamV3.h"
  42 #include "AliDAQ.h"
  43
  44 #include "AliEMCALRecParam.h"
  45 #include "AliEMCALLoader.h"
  46 #include "AliEMCALGeometry.h"
  47 class AliEMCALDigitizer;
  48 #include "AliEMCALDigit.h"
  49 #include "AliEMCAL.h"
  50 #include "AliCaloCalibPedestal.h"
  51
  52 ClassImp(AliEMCALRawUtils)
  53
  54 // Signal shape parameters
  55 Int_t    AliEMCALRawUtils::fgTimeBins = 256; // number of sampling bins of the raw RO signal (we typically use 15-50; theoretical max is 1k+)
  56 Double_t AliEMCALRawUtils::fgTimeBinWidth  = 100E-9 ; // each sample is 100 ns
  57 Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ;   // 15 time bins ~ 1.5 musec
  58
  59 // some digitization constants
  60 Int_t    AliEMCALRawUtils::fgThreshold = 1;
  61 Int_t    AliEMCALRawUtils::fgDDLPerSuperModule = 2;  // 2 ddls per SuperModule
  62 Int_t    AliEMCALRawUtils::fgPedestalValue = 32;     // pedestal value for digits2raw
  63 Double_t AliEMCALRawUtils::fgFEENoise = 3.;          // 3 ADC channels of noise (sampled)
  64
  65 AliEMCALRawUtils::AliEMCALRawUtils()
  66   : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
  67     fNPedSamples(0), fGeom(0), fOption("")
  68 {
  69
  70   //These are default parameters.
  71   //Can be re-set from without with setter functions
  72   fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits)
  73   fOrder = 2;                         // order of gamma fn
  74   fTau = 2.35;                        // in units of timebin, from CERN 2007 testbeam
  75   fNoiseThreshold = 3; // 3 ADC counts is approx. noise level
  76   fNPedSamples = 4;    // less than this value => likely pedestal samples
  77
  78   //Get Mapping RCU files from the AliEMCALRecParam
  79   const TObjArray* maps = AliEMCALRecParam::GetMappings();
  80   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
  81
  82   for(Int_t i = 0; i < 4; i++) {
  83     fMapping[i] = (AliAltroMapping*)maps->At(i);
  84   }
  85
  86   //To make sure we match with the geometry in a simulation file,
  87   //let's try to get it first.  If not, take the default geometry
  88   AliRunLoader *rl = AliRunLoader::Instance();
  89   if(!rl) AliError("Cannot find RunLoader!");
  90   if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) {
  91     fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
  92   } else {
  93     AliInfo(Form("Using default geometry in raw reco"));
  94     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
  95   }
  96
  97   if(!fGeom) AliFatal(Form("Could not get geometry!"));
  98
  99 }
 100
 101 //____________________________________________________________________________
 102 AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry)
 103   : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
 104     fNPedSamples(0), fGeom(pGeometry), fOption("")
 105 {
 106   //
 107   // Initialize with the given geometry - constructor required by HLT
 108   // HLT does not use/support AliRunLoader(s) instances
 109   // This is a minimum intervention solution
 110   // Comment by MPloskon@lbl.gov
 111   //
 112
 113   //These are default parameters.
 114   //Can be re-set from without with setter functions
 115   fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits)
 116   fOrder = 2;                         // order of gamma fn
 117   fTau = 2.35;                        // in units of timebin, from CERN 2007 testbeam
 118   fNoiseThreshold = 3; // 3 ADC counts is approx. noise level
 119   fNPedSamples = 4;    // less than this value => likely pedestal samples
 120
 121   //Get Mapping RCU files from the AliEMCALRecParam
 122   const TObjArray* maps = AliEMCALRecParam::GetMappings();
 123   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
 124
 125   for(Int_t i = 0; i < 4; i++) {
 126     fMapping[i] = (AliAltroMapping*)maps->At(i);
 127   }
 128
 129   if(!fGeom) AliFatal(Form("Could not get geometry!"));
 130
 131 }
 132
 133 //____________________________________________________________________________
 134 AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU)
 135   : TObject(),
 136     fHighLowGainFactor(rawU.fHighLowGainFactor),
 137     fOrder(rawU.fOrder),
 138     fTau(rawU.fTau),
 139     fNoiseThreshold(rawU.fNoiseThreshold),
 140     fNPedSamples(rawU.fNPedSamples),
 141     fGeom(rawU.fGeom),
 142     fOption(rawU.fOption)
 143 {
 144   //copy ctor
 145   fMapping[0] = rawU.fMapping[0];
 146   fMapping[1] = rawU.fMapping[1];
 147   fMapping[2] = rawU.fMapping[2];
 148   fMapping[3] = rawU.fMapping[3];
 149 }
 150
 151 //____________________________________________________________________________
 152 AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU)
 153 {
 154   //assignment operator
 155
 156   if(this != &rawU) {
 157     fHighLowGainFactor = rawU.fHighLowGainFactor;
 158     fOrder = rawU.fOrder;
 159     fTau = rawU.fTau;
 160     fNoiseThreshold = rawU.fNoiseThreshold;
 161     fNPedSamples = rawU.fNPedSamples;
 162     fGeom = rawU.fGeom;
 163     fOption = rawU.fOption;
 164     fMapping[0] = rawU.fMapping[0];
 165     fMapping[1] = rawU.fMapping[1];
 166     fMapping[2] = rawU.fMapping[2];
 167     fMapping[3] = rawU.fMapping[3];
 168   }
 169
 170   return *this;
 171
 172 }
 173
 174 //____________________________________________________________________________
 175 AliEMCALRawUtils::~AliEMCALRawUtils() {
 176   //dtor
 177
 178 }
 179
 180 //____________________________________________________________________________
 181 void AliEMCALRawUtils::Digits2Raw()
 182 {
 183   // convert digits of the current event to raw data
 184
 185   AliRunLoader *rl = AliRunLoader::Instance();
 186   AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 187
 188   // get the digits
 189   loader->LoadDigits("EMCAL");
 190   loader->GetEvent();
 191   TClonesArray* digits = loader->Digits() ;
 192
 193   if (!digits) {
 194     Warning("Digits2Raw", "no digits found !");
 195     return;
 196   }
 197
 198   static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
 199   AliAltroBuffer* buffers[nDDL];
 200   for (Int_t i=0; i < nDDL; i++)
 201     buffers[i] = 0;
 202
 203   TArrayI adcValuesLow(fgTimeBins);
 204   TArrayI adcValuesHigh(fgTimeBins);
 205
 206   // loop over digits (assume ordered digits)
 207   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
 208     AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
 209     if (digit->GetAmp() < fgThreshold)
 210       continue;
 211
 212     //get cell indices
 213     Int_t nSM = 0;
 214     Int_t nIphi = 0;
 215     Int_t nIeta = 0;
 216     Int_t iphi = 0;
 217     Int_t ieta = 0;
 218     Int_t nModule = 0;
 219     fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
 220     fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
 221
 222     //Check which is the RCU, 0 or 1, of the cell.
 223     Int_t iRCU = -111;
 224     //RCU0
 225     if (0<=iphi&&iphi<8) iRCU=0; // first cable row
 226     else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half;
 227     //second cable row
 228     //RCU1
 229     else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half;
 230     //second cable row
 231     else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
 232
 233     if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
 234
 235     if (iRCU<0)
 236       Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU);
 237
 238     //Which DDL?
 239     Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU;
 240     if (iDDL >= nDDL)
 241       Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
 242
 243     if (buffers[iDDL] == 0) {
 244       // open new file and write dummy header
 245       TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
 246       //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C
 247       Int_t iRCUside=iRCU+(nSM%2)*2;
 248       //iRCU=0 and even (0) SM -> RCU0A.data   0
 249       //iRCU=1 and even (0) SM -> RCU1A.data   1
 250       //iRCU=0 and odd  (1) SM -> RCU0C.data   2
 251       //iRCU=1 and odd  (1) SM -> RCU1C.data   3
 252       //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl;
 253       buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]);
 254       buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
 255     }
 256
 257     // out of time range signal (?)
 258     if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
 259       AliInfo("Signal is out of time range.\n");
 260       buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp());
 261       buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() );  // time bin
 262       buffers[iDDL]->FillBuffer(3);          // bunch length
 263       buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM);  // trailer
 264       // calculate the time response function
 265     } else {
 266       Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ;
 267       if (lowgain)
 268         buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow.GetArray(), fgThreshold);
 269       else
 270         buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh.GetArray(), fgThreshold);
 271     }
 272   }
 273
 274   // write headers and close files
 275   for (Int_t i=0; i < nDDL; i++) {
 276     if (buffers[i]) {
 277       buffers[i]->Flush();
 278       buffers[i]->WriteDataHeader(kFALSE, kFALSE);
 279       delete buffers[i];
 280     }
 281   }
 282
 283   loader->UnloadDigits();
 284 }
 285
 286 //____________________________________________________________________________
 287 void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, AliCaloCalibPedestal* pedbadmap)
 288 {
 289   // convert raw data of the current event to digits
 290
 291   digitsArr->Clear();
 292
 293   if (!digitsArr) {
 294     Error("Raw2Digits", "no digits found !");
 295     return;
 296   }
 297   if (!reader) {
 298     Error("Raw2Digits", "no raw reader found !");
 299     return;
 300   }
 301
 302   AliCaloRawStreamV3 in(reader,"EMCAL",fMapping);
 303   // Select EMCAL DDL's;
 304   reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL
 305
 306   //Updated fitting routine from 2007 beam test takes into account
 307   //possibility of two peaks in data and selects first one for fitting
 308   //Also sets some of the starting parameters based on the shape of the
 309   //given raw signal being fit
 310
 311   TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
 312   signalF->SetParameters(10.,5.,fTau,fOrder,0.); //set all defaults once, just to be safe
 313   signalF->SetParNames("amp","t0","tau","N","ped");
 314   signalF->FixParameter(2,fTau); // tau in units of time bin
 315   signalF->FixParameter(3,fOrder); // order
 316
 317   Int_t id =  -1;
 318   Float_t time = 0. ;
 319   Float_t amp = 0. ;
 320   Float_t ped = 0. ;
 321   Float_t ampEstimate  = 0;
 322   Float_t timeEstimate = 0;
 323   Float_t pedEstimate = 0;
 324   Int_t i = 0;
 325   Int_t startBin = 0;
 326
 327   //Graph to hold data we will fit (should be converted to an array
 328   //later to speed up processing
 329   TGraph * gSig = new TGraph(GetRawFormatTimeBins());
 330
 331   Int_t lowGain = 0;
 332   Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref.
 333
 334   // start loop over input stream
 335   while (in.NextDDL()) {
 336     while (in.NextChannel()) {
 337
 338       //Check if the signal  is high or low gain and then do the fit,
 339       //if it  is from TRU do not fit
 340       caloFlag = in.GetCaloFlag();
 341       if (caloFlag != 0 && caloFlag != 1) continue;
 342
 343       //Do not fit bad channels
 344       if(pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) {
 345         //printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow());
 346         continue;
 347       }
 348
 349       // There can be zero-suppression in the raw data,
 350       // so set up the TGraph in advance
 351       for (i=0; i < GetRawFormatTimeBins(); i++) {
 352         gSig->SetPoint(i, i , -1); // init to out-of-range values
 353       }
 354
 355       Int_t maxTimeBin = 0;
 356       Int_t min = 0x3ff; // init to 10-bit max
 357       Int_t max = 0; // init to 10-bit min
 358       while (in.NextBunch()) {
 359
 360         const UShort_t *sig = in.GetSignals();
 361         startBin = in.GetStartTimeBin();
 362         if (maxTimeBin < startBin) {
 363           maxTimeBin = startBin; // timebins come in reverse order
 364         }
 365         if (maxTimeBin < 0 || maxTimeBin >= GetRawFormatTimeBins()) {
 366           AliWarning(Form("Invalid time bin %d",maxTimeBin));
 367           maxTimeBin = GetRawFormatTimeBins();
 368         }
 369
 370         for (i = 0; i < in.GetBunchLength(); i++) {
 371           time = startBin--;
 372           gSig->SetPoint((Int_t)time, time, (Double_t) sig[i]) ;
 373           if (max < sig[i]) max = sig[i];
 374           if (min > sig[i]) min = sig[i];
 375
 376         }
 377       } // loop over bunches
 378
 379       gSig->Set(maxTimeBin+1); // set actual max size of TGraph
 380
 381       //Initialize the variables, do not keep previous values.
 382       // not really necessary to reset all of them (only amp and time at the moment), but better safe than sorry
 383       amp  = -1 ;
 384       time = -1 ;
 385       ped = -1;
 386       ampEstimate  = -1 ;
 387       timeEstimate = -1 ;
 388       pedEstimate = -1;
 389       if ( (max - min) > fNoiseThreshold) {
 390         FitRaw(gSig, signalF, maxTimeBin, amp, time, ped,
 391                ampEstimate, timeEstimate, pedEstimate);
 392       }
 393
 394       if ( amp>0 && amp<2000 && time>0 && time<(maxTimeBin*GetRawFormatTimeBinWidth()) ) {  //check both high and low end of amplitude result, and time
 395         //2000 is somewhat arbitrary - not nice with magic numbers in the code..
 396         id =  fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
 397         lowGain = in.IsLowGain();
 398
 399         // check fit results: should be consistent with initial estimates
 400         // more magic numbers, but very loose cuts, for now..
 401         // We have checked that amp an time values are positive so division for assymmetry
 402         // calculation should be OK/safe
 403         Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
 404         if ( (TMath::Abs(ampAsymm) > 0.1) ||
 405              (TMath::Abs(time - timeEstimate) > 2*GetRawFormatTimeBinWidth()) ) {
 406           AliDebug(2,Form("Fit results ped %f amp %f time %f not consistent with expectations ped %f max-ped %f time %d",
 407                       ped, amp, time, pedEstimate, ampEstimate, timeEstimate));
 408
 409           // what should do we do then? skip this channel or assign the simple estimate?
 410           // for now just overwrite the fit results with the simple estimate
 411           amp = ampEstimate;
 412           time = timeEstimate;
 413         }
 414
 415         AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
 416         // printf("Added tower: SM %d, row %d, column %d, amp %3.2f\n",in.GetModule(), in.GetRow(), in.GetColumn(),amp);
 417         // round off amplitude value to nearest integer
 418         AddDigit(digitsArr, id, lowGain, TMath::Nint(amp), time);
 419       }
 420
 421       // Reset graph
 422       for (Int_t index = 0; index < gSig->GetN(); index++) {
 423         gSig->SetPoint(index, index, -1) ;
 424       }
 425       // Reset starting parameters for fit function
 426       signalF->SetParameters(10.,5.,fTau,fOrder,0.); //reset all defaults just to be safe
 427
 428    } // end while over channel
 429   } //end while over DDL's, of input stream
 430
 431   delete signalF ;
 432   delete gSig;
 433
 434   return ;
 435 }
 436
 437 //____________________________________________________________________________
 438 void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) {
 439   //
 440   // Add a new digit.
 441   // This routine checks whether a digit exists already for this tower
 442   // and then decides whether to use the high or low gain info
 443   //
 444   // Called by Raw2Digits
 445
 446   AliEMCALDigit *digit = 0, *tmpdigit = 0;
 447   TIter nextdigit(digitsArr);
 448   while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) {
 449     if (tmpdigit->GetId() == id)
 450       digit = tmpdigit;
 451   }
 452
 453   if (!digit) { // no digit existed for this tower; create one
 454     if (lowGain && amp > fgkOverflowCut)
 455       amp = Int_t(fHighLowGainFactor * amp);
 456     Int_t idigit = digitsArr->GetEntries();
 457     new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ;
 458   }
 459   else { // a digit already exists, check range
 460          // (use high gain if signal < cut value, otherwise low gain)
 461     if (lowGain) { // new digit is low gain
 462       if (digit->GetAmp() > fgkOverflowCut) {  // use if stored digit is out of range
 463         digit->SetAmp(Int_t(fHighLowGainFactor * amp));
 464         digit->SetTime(time);
 465       }
 466     }
 467     else if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range
 468       digit->SetAmp(amp);
 469       digit->SetTime(time);
 470     }
 471   }
 472 }
 473
 474 //____________________________________________________________________________
 475 void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & ped, Float_t & ampEstimate, Float_t & timeEstimate, Float_t & pedEstimate, const Float_t cut) const
 476 {
 477   // Fits the raw signal time distribution; from AliEMCALGetter
 478   // last argument: Float_t cut = 0.0; // indicating how much of amplitude w.r.t. max value fit should be above noise and pedestal
 479
 480   // initialize return values
 481   amp = 0;
 482   time = 0;
 483   ped = 0;
 484   ampEstimate = 0;
 485   timeEstimate = 0;
 486   pedEstimate = 0;
 487
 488   // 0th step: remove plateau / overflow candidates
 489   // before trying to estimate amplitude, search for maxima etc.
 490   //
 491   Int_t nOrig = gSig->GetN(); // number of samples before we remove any overflows
 492   // Values for readback from input graph
 493   Double_t ttime = 0;
 494   Double_t signal = 0;
 495
 496   /*
 497   // start: tmp dump of all values
 498   for (Int_t i=0; i<gSig->GetN(); i++) {
 499     gSig->GetPoint(i, ttime, signal) ; // get values
 500     printf("orig: i %d, time %f, signal %f\n",i, ttime, signal);
 501   }
 502   // end: tmp dump of all values
 503   */
 504
 505   // start from back of TGraph since RemovePoint will downshift indices
 506   for (Int_t i=nOrig-1; i>=0; i--) {
 507     gSig->GetPoint(i, ttime, signal) ; // get values
 508     if (signal >= (pedEstimate + fgkOverflowCut) ) {
 509       gSig->RemovePoint(i);
 510     }
 511   }
 512
 513   // 1st step: we try to estimate the pedestal value
 514   Int_t nPed = 0;
 515   for (Int_t index = 0; index < gSig->GetN(); index++) {
 516     gSig->GetPoint(index, ttime, signal) ;
 517     // ttime < fNPedsamples used for pedestal estimate;
 518     // ttime >= fNPedSamples used for signal checks
 519     if (signal >= 0 && ttime<fNPedSamples) { // valid value
 520       pedEstimate += signal;
 521       nPed++;
 522     }
 523   }
 524
 525   if (nPed > 0)
 526     pedEstimate /= nPed;
 527   else {
 528     //AliWarning("Could not determine pedestal");
 529     AliDebug(1,"Could not determine pedestal");
 530     pedEstimate = 0; // good estimate for ZeroSupp data (non ZS data should have no problem with pedestal estimate)
 531   }
 532
 533   // 2nd step: we look through the rest of the time-bins/ADC values and
 534   // see if we have something that looks like a signal.
 535   // We look for a first local maxima, as well as for a global maxima
 536   Int_t locMaxFound = 0;
 537   Int_t locMaxId = 0; // time-bin index at first local max
 538   Float_t locMaxSig = -1; // actual local max value
 539   Int_t globMaxId = 0; // time-bin index at global max
 540   Float_t globMaxSig = -1; // actual global max value
 541   // We will also look for any values that look like they are in overflow region
 542   for (Int_t i=0; i<gSig->GetN(); i++) {
 543     gSig->GetPoint(i, ttime, signal) ; // get values
 544
 545     // ttime < fNPedsamples used for pedestal estimate;
 546     // ttime >= fNPedSamples used for signal checks
 547     if (ttime >= fNPedSamples) {
 548
 549       // look for first local maximum signal=ADC value
 550       if (!locMaxFound && signal > locMaxSig) {
 551         locMaxId = i;
 552         locMaxSig = signal;
 553       }
 554       else if ( locMaxSig > (pedEstimate + fNoiseThreshold) ) {
 555         // we enter this condition after signal<=max, but previous
 556         // max value was large enough. I.e. at least a significant local
 557         // maxima has been found (just before)
 558         locMaxFound = 1;
 559       }
 560
 561       // also check for global maximum..
 562       if (signal > globMaxSig) {
 563         globMaxId = i;
 564         globMaxSig = signal;
 565       }
 566     } // ttime check
 567   } // end for-loop over samples after pedestal
 568
 569   // OK, we have looked through the signal spectra, let's see if we should try to make the fit
 570   ampEstimate = locMaxSig - pedEstimate; // estimate using first local maxima
 571   if ( ampEstimate > fNoiseThreshold ) { // else it's just noise
 572
 573     //Check that the local maximum we will use is not at the end or beginning of time sample range
 574     Double_t timeMax = -1;
 575     Int_t iMax = locMaxId;
 576     gSig->GetPoint(locMaxId, timeMax, signal) ;
 577     if (timeMax < 2 || timeMax > lastTimeBin-1) { // lastTimeBin is the lowest kept time-sample; current (Dec 2009) case
 578       //    if (timeMax < 2 || timeMax > lastTimeBin-2) { // for when lastTimeBin is the lowest read-out time-sample, future (2010) case
 579       AliDebug(1,Form("Skip fit, maximum of the sample close to the edges : timeMax %3.2f, ampEstimate %3.2f",timeMax, ampEstimate));
 580       return;
 581     }
 582
 583     // Check if the local and global maximum disagree
 584     if (locMaxId != globMaxId) {
 585       AliDebug(1,Form("Warning, local first maximum %d does not agree with global maximum %d\n", locMaxId, globMaxId));
 586       return;
 587     }
 588
 589     // Get the maximum and find the lowest timebin (tailmin) where the ADC value is not
 590     // significantly different from the pedestal
 591     // first lower times edge a.k.a. tailmin
 592     Int_t tailMin = 0;
 593     Double_t tmptime = 0;
 594     for (Int_t i=iMax-1; i > 0; i--) {
 595       gSig->GetPoint(i, tmptime, signal) ;
 596       if((signal-pedEstimate) < fNoiseThreshold){
 597         tailMin = i;
 598         break;
 599       }
 600     }
 601     // then same exercise for the higher times edge a.k.a. tailmax
 602     Int_t tailMax = lastTimeBin;
 603     for (Int_t i=iMax+1; i < gSig->GetN(); i++) {
 604       gSig->GetPoint(i, tmptime, signal) ;
 605       if ((signal-pedEstimate) <= (ampEstimate*cut + fNoiseThreshold)) { // stop fit at cut-fraction of amplitude above noise-threshold (cut>0 would mean avoid the pulse shape falling edge)
 606         tailMax = i;
 607         break;
 608       }
 609     }
 610
 611     // remove all points which are not in the distribution around maximum
 612     // i.e. up to tailmin, and from tailmax
 613     if ( tailMax != (gSig->GetN()-1) ){ // else nothing to remove
 614       nOrig = gSig->GetN(); // can't use GetN call in for loop below since gSig size changes..
 615       for(int j = tailMax; j < nOrig; j++) gSig->RemovePoint(tailMax);
 616     }
 617     for(int j = 0; j<=tailMin; j++) gSig->RemovePoint(0);
 618
 619     if(gSig->GetN() < 3) {
 620       AliDebug(2,Form("Skip fit, number of entries in sample smaller than number of fitting parameters: in sample %d, fitting param 3",
 621                       gSig->GetN() ));
 622       return;
 623     }
 624
 625     timeEstimate = timeMax * GetRawFormatTimeBinWidth();
 626
 627     // determine what the valid fit range is
 628     Double_t minFit = 9999;
 629     Double_t maxFit = 0;
 630     for (Int_t i=0; i < gSig->GetN(); i++) {
 631       gSig->GetPoint(i, ttime, signal);
 632       if (minFit > ttime) minFit=ttime;
 633       if (maxFit < ttime) maxFit=ttime;
 634       //debug: printf("no tail: i %d, time %f, signal %f\n",i, ttime, signal);
 635     }
 636     signalF->SetRange(minFit, maxFit);
 637
 638     signalF->FixParameter(4, pedEstimate) ;
 639     signalF->SetParameter(1, timeMax);
 640     signalF->SetParameter(0, ampEstimate);
 641
 642     gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
 643
 644     // assign fit results
 645     amp = signalF->GetParameter(0);
 646     time = signalF->GetParameter(1) * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger?
 647     ped = signalF->GetParameter(4);
 648
 649     //BEG YS alternative methods to calculate the amplitude
 650     Double_t * ymx = gSig->GetX() ;
 651     Double_t * ymy = gSig->GetY() ;
 652     const Int_t kN = 3 ;
 653     Double_t ymMaxX[kN] = {0., 0., 0.} ;
 654     Double_t ymMaxY[kN] = {0., 0., 0.} ;
 655     Double_t ymax = 0. ;
 656       // find the maximum amplitude
 657     Int_t ymiMax = 0 ;
 658     for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
 659       if (ymy[ymi] > ymMaxY[0] ) {
 660         ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
 661         ymMaxX[0] = ymx[ymi] ;
 662         ymiMax = ymi ;
 663       }
 664     }
 665       // find the maximum by fitting a parabola through the max and the two adjacent samples
 666     if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
 667       ymMaxY[1] = ymy[ymiMax+1] ;
 668       ymMaxY[2] = ymy[ymiMax-1] ;
 669       ymMaxX[1] = ymx[ymiMax+1] ;
 670       ymMaxX[2] = ymx[ymiMax-1] ;
 671       if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
 672           //fit a parabola through the 3 points y= a+bx+x*x*x
 673         Double_t sy = 0 ;
 674         Double_t sx = 0 ;
 675         Double_t sx2 = 0 ;
 676         Double_t sx3 = 0 ;
 677         Double_t sx4 = 0 ;
 678         Double_t sxy = 0 ;
 679         Double_t sx2y = 0 ;
 680         for (Int_t i = 0; i < kN ; i++) {
 681           sy += ymMaxY[i] ;
 682           sx += ymMaxX[i] ;
 683           sx2 += ymMaxX[i]*ymMaxX[i] ;
 684           sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
 685           sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
 686           sxy += ymMaxX[i]*ymMaxY[i] ;
 687           sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ;
 688         }
 689         Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2);
 690         Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
 691         Double_t c  = cN / cD ;
 692         Double_t b  = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
 693         Double_t a  = (sy-b*sx-c*sx2)/kN  ;
 694         Double_t xmax = -b/(2*c) ;
 695         ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
 696       }
 697     }
 698
 699     Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ;
 700     if (diff > 0.1)
 701       amp = ymMaxY[0] ;
 702
 703       //END YS
 704
 705   } // ampEstimate > fNoiseThreshold
 706   return;
 707 }
 708 //__________________________________________________________________
 709 Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
 710 {
 711   // Matches version used in 2007 beam test
 712   //
 713   // Shape of the electronics raw reponse:
 714   // It is a semi-gaussian, 2nd order Gamma function of the general form
 715   //
 716   // xx = (t - t0 + tau) / tau  [xx is just a convenient help variable]
 717   // F = A * (xx**N * exp( N * ( 1 - xx) )   for xx >= 0
 718   // F = 0                                   for xx < 0
 719   //
 720   // parameters:
 721   // A:   par[0]   // Amplitude = peak value
 722   // t0:  par[1]
 723   // tau: par[2]
 724   // N:   par[3]
 725   // ped: par[4]
 726   //
 727   Double_t signal ;
 728   Double_t tau =par[2];
 729   Double_t n =par[3];
 730   Double_t ped = par[4];
 731   Double_t xx = ( x[0] - par[1] + tau ) / tau ;
 732
 733   if (xx <= 0)
 734     signal = ped ;
 735   else {
 736     signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
 737   }
 738   return signal ;
 739 }
 740
 741 //__________________________________________________________________
 742 Bool_t AliEMCALRawUtils::RawSampledResponse(
 743 const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const
 744 {
 745   // for a start time dtime and an amplitude damp given by digit,
 746   // calculates the raw sampled response AliEMCAL::RawResponseFunction
 747
 748   Bool_t lowGain = kFALSE ;
 749
 750   // A:   par[0]   // Amplitude = peak value
 751   // t0:  par[1]
 752   // tau: par[2]
 753   // N:   par[3]
 754   // ped: par[4]
 755
 756   TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
 757   signalF.SetParameter(0, damp) ;
 758   signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ;
 759   signalF.SetParameter(2, fTau) ;
 760   signalF.SetParameter(3, fOrder);
 761   signalF.SetParameter(4, fgPedestalValue);
 762
 763   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 764     Double_t signal = signalF.Eval(iTime) ;
 765
 766     // Next lines commeted for the moment but in principle it is not necessary to add
 767     // extra noise since noise already added at the digits level.
 768
 769     //According to Terry Awes, 13-Apr-2008
 770     //add gaussian noise in quadrature to each sample
 771     //Double_t noise = gRandom->Gaus(0.,fgFEENoise);
 772     //signal = sqrt(signal*signal + noise*noise);
 773
 774     // March 17,09 for fast fit simulations by Alexei Pavlinov.
 775     // Get from PHOS analysis. In some sense it is open questions.
 776     //Double_t noise = gRandom->Gaus(0.,fgFEENoise);
 777     //signal += noise;
 778
 779     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
 780     if ( adcH[iTime] > fgkRawSignalOverflow ){  // larger than 10 bits
 781       adcH[iTime] = fgkRawSignalOverflow ;
 782       lowGain = kTRUE ;
 783     }
 784
 785     signal /= fHighLowGainFactor;
 786
 787     adcL[iTime] =  static_cast<Int_t>(signal + 0.5) ;
 788     if ( adcL[iTime] > fgkRawSignalOverflow)  // larger than 10 bits
 789       adcL[iTime] = fgkRawSignalOverflow ;
 790   }
 791   return lowGain ;
 792 }