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 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.8  2007/12/05 02:30:51  jklay
  21  * modification to read Altro mappings into AliEMCALRecParam and pass to AliEMCALRawUtils from AliEMCALReconstructor; add option to AliEMCALRawUtils to set old RCU format (for testbeam) or not
  22  *
  23  * Revision 1.7  2007/11/14 15:51:46  gustavo
  24  * Take out few unnecessary prints
  25  *
  26  * Revision 1.6  2007/11/01 01:23:51  mvl
  27  * Removed call to SetOldRCUFormat, which is only needed for testbeam data
  28  *
  29  * Revision 1.5  2007/11/01 01:20:33  mvl
  30  * Further improvement of peak finding; more robust fit
  31  *
  32  * Revision 1.4  2007/10/31 17:15:24  mvl
  33  * Fixed bug in raw data unpacking; Added pedestal to signal fit; Added logic to deal with high/low gain
  34  *
  35  * Revision 1.3  2007/09/27 08:36:46  mvl
  36  * More robust setting of fit range in FitRawSignal (P. Hristov)
  37  *
  38  * Revision 1.2  2007/09/03 20:55:35  jklay
  39  * EMCAL e-by-e reconstruction methods from Cvetan
  40  *
  41  * Revision 1.1  2007/03/17 19:56:38  mvl
  42  * Moved signal shape routines from AliEMCAL to separate class AliEMCALRawUtils to streamline raw data reconstruction code.
  43  * */
  44
  45 //*-- Author: Marco van Leeuwen (LBL)
  46 #include "AliEMCALRawUtils.h"
  47
  48 #include "TF1.h"
  49 #include "TGraph.h"
  50 #include "TSystem.h"
  51
  52 #include "AliLog.h"
  53 #include "AliRunLoader.h"
  54 #include "AliCaloAltroMapping.h"
  55 #include "AliAltroBuffer.h"
  56 #include "AliRawReader.h"
  57 #include "AliCaloRawStream.h"
  58 #include "AliDAQ.h"
  59
  60 #include "AliEMCALRecParam.h"
  61 #include "AliEMCALLoader.h"
  62 #include "AliEMCALGeometry.h"
  63 #include "AliEMCALDigitizer.h"
  64 #include "AliEMCALDigit.h"
  65
  66
  67 ClassImp(AliEMCALRawUtils)
  68
  69 // Signal shape parameters
  70 Int_t    AliEMCALRawUtils::fgOrder         = 2 ;      // Order of gamma function
  71 Double_t AliEMCALRawUtils::fgTimeBinWidth  = 100E-9 ; // each sample is 100 ns
  72 Double_t AliEMCALRawUtils::fgTau         = 235E-9 ;   // 235 ns (from CERN testbeam; not very accurate)
  73 Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ;   // 15 time bins ~ 1.5 musec
  74
  75 // some digitization constants
  76 Int_t    AliEMCALRawUtils::fgThreshold = 1;
  77 Int_t    AliEMCALRawUtils::fgDDLPerSuperModule = 2;  // 2 ddls per SuperModule
  78
  79 AliEMCALRawUtils::AliEMCALRawUtils()
  80   : fHighLowGainFactor(0.), fOption("")
  81 {
  82   fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits)
  83 }
  84 //____________________________________________________________________________
  85 AliEMCALRawUtils::~AliEMCALRawUtils() {
  86 }
  87 //____________________________________________________________________________
  88 void AliEMCALRawUtils::Digits2Raw(AliAltroMapping **mapping)
  89 {
  90   // convert digits of the current event to raw data
  91
  92   AliRunLoader *rl = AliRunLoader::GetRunLoader();
  93   AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
  94
  95   // get the digits
  96   loader->LoadDigits("EMCAL");
  97   loader->GetEvent();
  98   TClonesArray* digits = loader->Digits() ;
  99
 100   if (!digits) {
 101     Warning("Digits2Raw", "no digits found !");
 102     return;
 103   }
 104
 105   // get the geometry
 106   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 107   if (!geom) {
 108     AliError(Form("No geometry found !"));
 109     return;
 110   }
 111
 112   static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
 113   AliAltroBuffer* buffers[nDDL];
 114   for (Int_t i=0; i < nDDL; i++)
 115     buffers[i] = 0;
 116
 117   Int_t adcValuesLow[fgkTimeBins];
 118   Int_t adcValuesHigh[fgkTimeBins];
 119
 120   // loop over digits (assume ordered digits)
 121   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
 122     AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
 123     if (digit->GetAmp() < fgThreshold)
 124       continue;
 125
 126     //get cell indices
 127     Int_t nSM = 0;
 128     Int_t nIphi = 0;
 129     Int_t nIeta = 0;
 130     Int_t iphi = 0;
 131     Int_t ieta = 0;
 132     Int_t nModule = 0;
 133     geom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
 134     geom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
 135
 136     //Check which is the RCU of the cell.
 137     Int_t iRCU = -111;
 138     //RCU0
 139     if (0<=iphi&&iphi<8) iRCU=0; // first cable row
 140     else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half;
 141     //second cable row
 142     //RCU1
 143     else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half;
 144     //second cable row
 145     else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
 146
 147     //Which DDL?
 148     Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU;
 149     if (iDDL >= nDDL)
 150       Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
 151
 152     if (buffers[iDDL] == 0) {
 153       // open new file and write dummy header
 154       TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
 155       buffers[iDDL] = new AliAltroBuffer(fileName.Data(),mapping[iRCU]);
 156       buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
 157     }
 158
 159     // out of time range signal (?)
 160     if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
 161       AliInfo("Signal is out of time range.\n");
 162       buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp());
 163       buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() );  // time bin
 164       buffers[iDDL]->FillBuffer(3);          // bunch length
 165       buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM);  // trailer
 166       // calculate the time response function
 167     } else {
 168       Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ;
 169       if (lowgain)
 170         buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold);
 171       else
 172         buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold);
 173     }
 174   }
 175
 176   // write headers and close files
 177   for (Int_t i=0; i < nDDL; i++) {
 178     if (buffers[i]) {
 179       buffers[i]->Flush();
 180       buffers[i]->WriteDataHeader(kFALSE, kFALSE);
 181       delete buffers[i];
 182     }
 183   }
 184   mapping[0]->Delete();
 185   mapping[1]->Delete();
 186   loader->UnloadDigits();
 187 }
 188
 189 //____________________________________________________________________________
 190 void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
 191                                   AliAltroMapping **mapping)
 192 {
 193   // convert raw data of the current event to digits
 194   AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 195   if (!geom) {
 196     AliError(Form("No geometry found !"));
 197     return;
 198   }
 199
 200   digitsArr->Clear();
 201
 202   if (!digitsArr) {
 203     Error("Raw2Digits", "no digits found !");
 204     return;
 205   }
 206   if (!reader) {
 207     Error("Raw2Digits", "no raw reader found !");
 208     return;
 209   }
 210
 211   AliCaloRawStream in(reader,"EMCAL",mapping);
 212   // Select EMCAL DDL's;
 213   reader->Select("EMCAL");
 214
 215   TString option = GetOption();
 216   if (option.Contains("OldRCUFormat"))
 217     in.SetOldRCUFormat(kTRUE); // Needed for testbeam data
 218   else
 219     in.SetOldRCUFormat(kFALSE);
 220
 221
 222   //Updated fitting routine from 2007 beam test takes into account
 223   //possibility of two peaks in data and selects first one for fitting
 224   //Also sets some of the starting parameters based on the shape of the
 225   //given raw signal being fit
 226
 227   TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
 228   signalF->SetParNames("amp","t0","tau","N","ped");
 229   signalF->SetParameter(2,2.35); // tau in units of time bin
 230   signalF->SetParLimits(2,2,-1);
 231   signalF->SetParameter(3,2); // order
 232   signalF->SetParLimits(3,2,-1);
 233
 234   Int_t id =  -1;
 235   Float_t time = 0. ;
 236   Float_t amp = 0. ;
 237
 238   //Graph to hold data we will fit (should be converted to an array
 239   //later to speed up processing
 240   TGraph * gSig = new TGraph(GetRawFormatTimeBins());
 241
 242   Int_t readOk = 1;
 243   Int_t lowGain = 0;
 244
 245   while (readOk && in.GetModule() < 0)
 246     readOk = in.Next();  // Go to first digit
 247
 248   Int_t col = 0;
 249   Int_t row = 0;
 250
 251   while (readOk) {
 252     id =  geom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
 253     lowGain = in.IsLowGain();
 254     Int_t maxTime = in.GetTime();  // timebins come in reverse order
 255     if (maxTime < 0 || maxTime >= GetRawFormatTimeBins()) {
 256       AliWarning(Form("Invalid time bin %d",maxTime));
 257       maxTime = GetRawFormatTimeBins();
 258     }
 259     gSig->Set(maxTime+1);
 260     // There is some kind of zero-suppression in the raw data,
 261     // so set up the TGraph in advance
 262     for (Int_t i=0; i < maxTime; i++) {
 263       gSig->SetPoint(i, i , 0);
 264     }
 265
 266     Int_t iTime = 0;
 267     do {
 268       if (in.GetTime() >= gSig->GetN()) {
 269           AliWarning("Too many time bins");
 270           gSig->Set(in.GetTime());
 271       }
 272       col = in.GetColumn();
 273       row = in.GetRow();
 274
 275       gSig->SetPoint(in.GetTime(), in.GetTime(), in.GetSignal()) ;
 276       if (in.GetTime() > maxTime)
 277         maxTime = in.GetTime();
 278       iTime++;
 279     } while ((readOk = in.Next()) && !in.IsNewHWAddress());
 280
 281     FitRaw(gSig, signalF, amp, time) ;
 282
 283     if (amp > 0) {
 284       AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
 285       //cout << "col " << col-40 << " row " << row-8 << " lowGain " << lowGain << " amp " << amp << endl;
 286       AddDigit(digitsArr, id, lowGain, (Int_t)amp, time);
 287     }
 288
 289     // Reset graph
 290     for (Int_t index = 0; index < gSig->GetN(); index++) {
 291       gSig->SetPoint(index, index, 0) ;
 292     }
 293   }; // EMCAL entries loop
 294
 295   delete signalF ;
 296   delete gSig;
 297
 298   return ;
 299 }
 300
 301 //____________________________________________________________________________
 302 void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) {
 303   //
 304   // Add a new digit.
 305   // This routine checks whether a digit exists already for this tower
 306   // and then decides whether to use the high or low gain info
 307   //
 308   // Called by Raw2Digits
 309
 310   AliEMCALDigit *digit = 0, *tmpdigit = 0;
 311
 312   TIter nextdigit(digitsArr);
 313   while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) {
 314     if (tmpdigit->GetId() == id)
 315       digit = tmpdigit;
 316   }
 317
 318   if (!digit) { // no digit existed for this tower; create one
 319     if (lowGain)
 320       amp = Int_t(fHighLowGainFactor * amp);
 321     Int_t idigit = digitsArr->GetEntries();
 322     new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ;
 323   }
 324   else { // a digit already exists, check range
 325          // (use high gain if signal < 800, otherwise low gain)
 326     if (lowGain) { // new digit is low gain
 327       if (digit->GetAmp() > 800) {  // use if stored digit is out of range
 328         digit->SetAmp(Int_t(fHighLowGainFactor * amp));
 329         digit->SetTime(time);
 330       }
 331     }
 332     else if (amp < 800) { // new digit is high gain; use if not out of range
 333       digit->SetAmp(amp);
 334       digit->SetTime(time);
 335     }
 336   }
 337 }
 338
 339 //____________________________________________________________________________
 340 void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_t & time)
 341 {
 342   // Fits the raw signal time distribution; from AliEMCALGetter
 343
 344   const Int_t kNoiseThreshold = 5;
 345   const Int_t kNPedSamples = 5;
 346   amp = time = 0. ;
 347   Double_t ped = 0;
 348   Int_t nPed = 0;
 349
 350   for (Int_t index = 0; index < kNPedSamples; index++) {
 351     Double_t ttime, signal;
 352     gSig->GetPoint(index, ttime, signal) ;
 353     if (signal > 0) {
 354       ped += signal;
 355       nPed++;
 356     }
 357   }
 358
 359   if (nPed > 0)
 360     ped /= nPed;
 361   else {
 362     AliWarning("Could not determine pedestal");
 363     ped = 10; // put some small value as first guess
 364   }
 365
 366   Int_t max_found = 0;
 367   Int_t i_max = 0;
 368   Float_t max = -1;
 369   Float_t max_fit = gSig->GetN();
 370   Float_t min_after_sig = 9999;
 371   Int_t tmin_after_sig = gSig->GetN();
 372   Int_t n_ped_after_sig = 0;
 373
 374   for (Int_t i=kNPedSamples; i < gSig->GetN(); i++) {
 375     Double_t ttime, signal;
 376     gSig->GetPoint(i, ttime, signal) ;
 377     if (!max_found && signal > max) {
 378       i_max = i;
 379       max = signal;
 380     }
 381     else if ( max > ped + kNoiseThreshold ) {
 382       max_found = 1;
 383       min_after_sig = signal;
 384       tmin_after_sig = i;
 385     }
 386     if (max_found) {
 387       if ( signal < min_after_sig) {
 388         min_after_sig = signal;
 389         tmin_after_sig = i;
 390       }
 391       if (i > tmin_after_sig + 5) {  // Two close peaks; end fit at minimum
 392         max_fit = tmin_after_sig;
 393         break;
 394       }
 395       if ( signal < ped + kNoiseThreshold)
 396         n_ped_after_sig++;
 397       if (n_ped_after_sig >= 5) {  // include 5 pedestal bins after peak
 398         max_fit = i;
 399         break;
 400       }
 401     }
 402   }
 403
 404   if ( max - ped > kNoiseThreshold ) { // else its noise
 405     AliDebug(2,Form("Fitting max %d ped %d", max, ped));
 406     signalF->SetRange(0,max_fit);
 407
 408     if(max-ped > 50)
 409       signalF->SetParLimits(2,1,3);
 410
 411     signalF->SetParameter(4, ped) ;
 412     signalF->SetParameter(1, i_max);
 413     signalF->SetParameter(0, max);
 414
 415     gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
 416     amp = signalF->GetParameter(0);
 417     time = signalF->GetParameter(1)*GetRawFormatTimeBinWidth() - fgTimeTrigger;
 418   }
 419   return;
 420 }
 421 //__________________________________________________________________
 422 Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
 423 {
 424   // Matches version used in 2007 beam test
 425   //
 426   // Shape of the electronics raw reponse:
 427   // It is a semi-gaussian, 2nd order Gamma function of the general form
 428   //
 429   // t' = (t - t0 + tau) / tau
 430   // F = A * t**N * exp( N * ( 1 - t) )   for t >= 0
 431   // F = 0                                for t < 0
 432   //
 433   // parameters:
 434   // A:   par[0]   // Amplitude = peak value
 435   // t0:  par[1]
 436   // tau: par[2]
 437   // N:   par[3]
 438   // ped: par[4]
 439   //
 440   Double_t signal ;
 441   Double_t tau =par[2];
 442   Double_t N =par[3];
 443   Double_t ped = par[4];
 444   Double_t xx = ( x[0] - par[1] + tau ) / tau ;
 445
 446   if (xx <= 0)
 447     signal = ped ;
 448   else {
 449     signal = ped + par[0] * TMath::Power(xx , N) * TMath::Exp(N * (1 - xx )) ;
 450   }
 451   return signal ;
 452 }
 453
 454 //__________________________________________________________________
 455 Bool_t AliEMCALRawUtils::RawSampledResponse(
 456 const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const
 457 {
 458   // for a start time dtime and an amplitude damp given by digit,
 459   // calculates the raw sampled response AliEMCAL::RawResponseFunction
 460
 461   const Int_t kRawSignalOverflow = 0x3FF ;
 462   const Int_t pedVal = 32;
 463   Bool_t lowGain = kFALSE ;
 464
 465   TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
 466   signalF.SetParameter(0, pedVal) ;
 467   signalF.SetParameter(1, damp) ;
 468   signalF.SetParameter(2, dtime + fgTimeTrigger) ;
 469
 470   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 471     Double_t time = iTime * GetRawFormatTimeBinWidth() ;
 472     Double_t signal = signalF.Eval(time) ;
 473     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
 474     if ( adcH[iTime] > kRawSignalOverflow ){  // larger than 10 bits
 475       adcH[iTime] = kRawSignalOverflow ;
 476       lowGain = kTRUE ;
 477     }
 478
 479     signal /= fHighLowGainFactor;
 480
 481     adcL[iTime] =  static_cast<Int_t>(signal + 0.5) ;
 482     if ( adcL[iTime] > kRawSignalOverflow)  // larger than 10 bits
 483       adcL[iTime] = kRawSignalOverflow ;
 484   }
 485   return lowGain ;
 486 }