EMCAL/AliEMCALRawUtils.cxx

   1 // -*- mode: c++ -*-
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /* $Id$ */
  18
  19
  20 //_________________________________________________________________________
  21 //  Utility Class for handling Raw data
  22 //  Does all transitions from Digits to Raw and vice versa,
  23 //  for simu and reconstruction
  24 //
  25 //  Note: the current version is still simplified. Only
  26 //    one raw signal per digit is generated; either high-gain or low-gain
  27 //    Need to add concurrent high and low-gain info in the future
  28 //    No pedestal is added to the raw signal.
  29 //*-- Author: Marco van Leeuwen (LBL)
  30
  31
  32 #include "AliEMCALRawUtils.h"
  33 #include "TF1.h"
  34 #include "TGraph.h"
  35 #include <TRandom.h>
  36 #include "AliRun.h"
  37 #include "AliRunLoader.h"
  38 #include "AliAltroBuffer.h"
  39 #include "AliRawReader.h"
  40 #include "AliCaloRawStreamV3.h"
  41 #include "AliDAQ.h"
  42 #include "AliEMCALRecParam.h"
  43 #include "AliEMCALLoader.h"
  44 #include "AliEMCALGeometry.h"
  45 #include "AliEMCALDigit.h"
  46 #include "AliEMCALRawDigit.h"
  47 #include "AliEMCAL.h"
  48 #include "AliCaloCalibPedestal.h"
  49 #include "AliCaloBunchInfo.h"
  50 #include "AliCaloFitResults.h"
  51 #include "AliEMCALTriggerRawDigitMaker.h"
  52 #include "AliEMCALTriggerSTURawStream.h"
  53 #include "AliEMCALTriggerData.h"
  54 #include "AliCaloConstants.h"
  55 #include "AliCaloRawAnalyzer.h"
  56 #include "AliCaloRawAnalyzerFactory.h"
  57
  58 using namespace CALO;
  59 using namespace EMCAL;
  60
  61 Double_t AliEMCALRawUtils::fgTimeTrigger  = 600E-9 ;   // the time of the trigger as approximately seen in the data
  62 Int_t    AliEMCALRawUtils::fgThreshold         = 1;
  63 Int_t    AliEMCALRawUtils::fgPedestalValue     = 0;  // pedestal value for digits2raw, default generate ZS data
  64 Double_t AliEMCALRawUtils::fgFEENoise          = 3.; // 3 ADC channels of noise (sampled)
  65
  66 ClassImp(AliEMCALRawUtils)
  67
  68
  69 AliEMCALRawUtils::AliEMCALRawUtils( Algo::fitAlgorithm fitAlgo) : fNoiseThreshold(3),
  70                                                                   fNPedSamples(4),
  71                                                                   fGeom(0),
  72                                                                   fOption(""),
  73                                                                   fRemoveBadChannels(kFALSE),
  74                                                                   fFittingAlgorithm(0),
  75                                                                   fTimeMin(-1.),
  76                                                                   fTimeMax(1.),
  77                                                                   fUseFALTRO(kTRUE),
  78                                                                   fRawAnalyzer(0),
  79                                                                   fTriggerRawDigitMaker(0x0)
  80 {
  81   SetFittingAlgorithm(fitAlgo);
  82   // SetFittingAlgorithm(  Algo::kLMSOffline);
  83
  84   //Get Mapping RCU files from the AliEMCALRecParam
  85   const TObjArray* maps = AliEMCALRecParam::GetMappings();
  86   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
  87
  88   for(Int_t i = 0; i < 4; i++) {
  89     fMapping[i] = (AliAltroMapping*)maps->At(i);
  90   }
  91
  92   //To make sure we match with the geometry in a simulation file,
  93   //let's try to get it first.  If not, take the default geometry
  94   AliRunLoader *rl = AliRunLoader::Instance();
  95   if (rl && rl->GetAliRun()) {
  96     AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
  97     if(emcal)fGeom = emcal->GetGeometry();
  98     else {
  99       AliDebug(1, Form("Using default geometry in raw reco"));
 100       fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
 101     }
 102
 103   } else {
 104     AliDebug(1, Form("Using default geometry in raw reco"));
 105     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
 106   }
 107
 108   if(!fGeom) AliFatal(Form("Could not get geometry!"));
 109
 110   fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();
 111
 112 }
 113
 114
 115 //____________________________________________________________________________
 116 AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry, Algo::fitAlgorithm fitAlgo) : //fHighLowGainFactor(16.),
 117                                                                                               //fOrder(2),
 118   //  fTau(2.35),
 119   fNoiseThreshold(3),
 120   fNPedSamples(4),
 121   fGeom(pGeometry),
 122   fOption(""),
 123   fRemoveBadChannels(kFALSE),fFittingAlgorithm(0),
 124   fTimeMin(-1.),fTimeMax(1.),
 125   fUseFALTRO(kTRUE),fRawAnalyzer(0),
 126   fTriggerRawDigitMaker(0x0)
 127 {
 128
 129
 130  // Initialize with the given geometry - constructor required by HLT
 131   // HLT does not use/support AliRunLoader(s) instances
 132   // This is a minimum intervention solution
 133   // Comment by MPloskon@lbl.gov
 134   SetFittingAlgorithm(fitAlgo);
 135   // SetFittingAlgorithm(  Algo::kLMSOffline);
 136   //Get Mapping RCU files from the AliEMCALRecParam
 137   const TObjArray* maps = AliEMCALRecParam::GetMappings();
 138   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
 139
 140   for(Int_t i = 0; i < 4; i++)
 141     {
 142       fMapping[i] = (AliAltroMapping*)maps->At(i);
 143     }
 144
 145   if(!fGeom) AliFatal(Form("Could not get geometry!"));
 146   fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();
 147 }
 148
 149
 150 //____________________________________________________________________________
 151 AliEMCALRawUtils::~AliEMCALRawUtils()
 152 {
 153   //dtor
 154 }
 155
 156
 157 //____________________________________________________________________________
 158 void AliEMCALRawUtils::Digits2Raw()
 159 {
 160   // convert digits of the current event to raw data
 161   AliRunLoader *rl = AliRunLoader::Instance();
 162   AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 163
 164   // get the digits
 165   loader->LoadDigits("EMCAL");
 166   loader->GetEvent();
 167   TClonesArray* digits = loader->Digits() ;
 168
 169   if (!digits) {
 170     Warning("Digits2Raw", "no digits found !");
 171     return;
 172   }
 173
 174   static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
 175   AliAltroBuffer* buffers[nDDL];
 176   for (Int_t i=0; i < nDDL; i++)
 177     buffers[i] = 0;
 178
 179   TArrayI adcValuesLow( TIMEBINS );
 180   TArrayI adcValuesHigh( TIMEBINS );
 181
 182   // loop over digits (assume ordered digits)
 183   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++)
 184     {
 185       AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
 186       if(!digit)
 187         {
 188           AliFatal("NULL Digit");
 189         }
 190       else
 191         {
 192           if (digit->GetAmplitude() < fgThreshold)
 193             {
 194               continue;
 195             }
 196           //get cell indices
 197           Int_t nSM = 0;
 198           Int_t nIphi = 0;
 199           Int_t nIeta = 0;
 200           Int_t iphi = 0;
 201           Int_t ieta = 0;
 202           Int_t nModule = 0;
 203           fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
 204           fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
 205
 206           //Check which is the RCU, 0 or 1, of the cell.
 207           Int_t iRCU = -111;
 208           //RCU0
 209       if (0<=iphi&&iphi<8) iRCU=0; // first cable row
 210       else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half;
 211       //second cable row
 212       //RCU1
 213       else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half;
 214       //second cable row
 215       else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
 216
 217       if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
 218
 219       if (iRCU<0)
 220         Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU);
 221
 222       //Which DDL?
 223       Int_t iDDL = NRCUSPERMODULE*nSM + iRCU;
 224       if (iDDL < 0 || iDDL >= nDDL){
 225         Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
 226       }
 227       else{
 228         if (buffers[iDDL] == 0) {
 229           // open new file and write dummy header
 230           TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
 231           //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C
 232           Int_t iRCUside=iRCU+(nSM%2)*2;
 233           //iRCU=0 and even (0) SM -> RCU0A.data   0
 234           //iRCU=1 and even (0) SM -> RCU1A.data   1
 235           //iRCU=0 and odd  (1) SM -> RCU0C.data   2
 236           //iRCU=1 and odd  (1) SM -> RCU1C.data   3
 237           //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl;
 238           buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]);
 239           buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
 240         }
 241
 242         // out of time range signal (?)
 243         if (digit->GetTimeR() >  TIMEBINMAX  ) {
 244           AliInfo("Signal is out of time range.\n");
 245           buffers[iDDL]->FillBuffer((Int_t)digit->GetAmplitude());
 246           buffers[iDDL]->FillBuffer( TIMEBINS );  // time bin
 247           buffers[iDDL]->FillBuffer(3);          // bunch length
 248           buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM);  // trailer
 249           // calculate the time response function
 250         } else {
 251           Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ;
 252           if (lowgain)
 253             buffers[iDDL]->WriteChannel(ieta, iphi, 0, TIMEBINS, adcValuesLow.GetArray(), fgThreshold);
 254           else
 255             buffers[iDDL]->WriteChannel(ieta,iphi, 1, TIMEBINS, adcValuesHigh.GetArray(), fgThreshold);
 256         }
 257       }// iDDL under the limits
 258     }//digit exists
 259   }//Digit loop
 260
 261   // write headers and close files
 262   for (Int_t i=0; i < nDDL; i++) {
 263     if (buffers[i]) {
 264       buffers[i]->Flush();
 265       buffers[i]->WriteDataHeader(kFALSE, kFALSE);
 266       delete buffers[i];
 267     }
 268   }
 269
 270   loader->UnloadDigits();
 271 }
 272
 273
 274 //____________________________________________________________________________
 275 void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf)
 276 {
 277   AliEMCALDigit *digit = 0, *tmpdigit = 0;
 278   TIter nextdigit(digitsArr);
 279
 280   while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit()))
 281     {
 282       if (tmpdigit->GetId() == id) digit = tmpdigit;
 283     }
 284
 285   if (!digit) { // no digit existed for this tower; create one
 286     Int_t type = AliEMCALDigit::kHG; // use enum in AliEMCALDigit
 287     if (lowGain)
 288       {
 289         amp *= HGLGFACTOR;
 290         type = AliEMCALDigit::kLGnoHG;
 291       }
 292
 293     Int_t idigit = digitsArr->GetEntries();
 294     new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit, chi2, ndf);
 295     AliDebug(2,Form("Add digit Id %d for the first time, type %d", id, type));
 296   }//digit added first time
 297   else
 298     { // a digit already exists, check range
 299                 // (use high gain if signal < cut value, otherwise low gain)
 300       if (lowGain)
 301         { // new digit is low gain
 302           if (digit->GetAmplitude() >  OVERFLOWCUT )
 303             {  // use if previously stored (HG) digit is out of range
 304               digit->SetAmplitude( HGLGFACTOR * amp);
 305               digit->SetTime(time);
 306               digit->SetType(AliEMCALDigit::kLG);
 307               AliDebug(2,Form("Add LG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType()));
 308             }
 309         }//new low gain digit
 310       else { // new digit is high gain
 311
 312         if (amp <  OVERFLOWCUT  )
 313           { // new digit is high gain; use if not out of range
 314             digit->SetAmplitude(amp);
 315             digit->SetTime(time);
 316             digit->SetType(AliEMCALDigit::kHG);
 317             AliDebug(2,Form("Add HG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType()));
 318           }
 319         else
 320           { // HG out of range, just change flag value to show that HG did exist
 321             digit->SetType(AliEMCALDigit::kLG);
 322             AliDebug(2,Form("Change LG digit to HG, ID %d, type %d", digit->GetId(), digit->GetType()));
 323           }
 324       }//new high gain digit
 325     }//digit existed replace it
 326 }
 327
 328
 329 //____________________________________________________________________________
 330 void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap, TClonesArray *digitsTRG, AliEMCALTriggerData* trgData)
 331 {
 332
 333   if(digitsArr) digitsArr->Clear("C");
 334   if (!digitsArr) { Error("Raw2Digits", "no digits found !");return;}
 335   if (!reader) {Error("Raw2Digits", "no raw reader found !");return;}
 336   AliEMCALTriggerSTURawStream inSTU(reader);
 337   AliCaloRawStreamV3 in(reader,"EMCAL",fMapping);
 338   reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL
 339   fTriggerRawDigitMaker->Reset();
 340   fTriggerRawDigitMaker->SetIO(reader, in, inSTU, digitsTRG, trgData);
 341   fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later
 342
 343   Int_t lowGain  = 0;
 344   Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref.
 345
 346   while (in.NextDDL())
 347     {
 348       //  fprintf(fp," TP1\n");
 349       while (in.NextChannel())
 350         {
 351           //      fprintf(fp," TP2\n");
 352           caloFlag = in.GetCaloFlag();
 353           if (caloFlag > 2) continue; // Work with ALTRO and FALTRO
 354           if(caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow()))
 355             {
 356               continue;
 357             }
 358           vector<AliCaloBunchInfo> bunchlist;
 359           while (in.NextBunch())
 360             {
 361               bunchlist.push_back( AliCaloBunchInfo(in.GetStartTimeBin(), in.GetBunchLength(), in.GetSignals() ) );
 362             }
 363           if (bunchlist.size() == 0) continue;
 364           if ( caloFlag < 2 )
 365             { // ALTRO
 366               Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
 367               lowGain  = in.IsLowGain();
 368               fRawAnalyzer->SetL1Phase( in.GetL1Phase() );
 369               AliCaloFitResults res =  fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2());
 370               if(res.GetAmp() >= fNoiseThreshold )
 371                 {
 372                   AddDigit(digitsArr, id, lowGain, res.GetAmp(),  res.GetTime(), res.GetChi2(),  res.GetNdf() );
 373                 }
 374             }//ALTRO
 375           else if(fUseFALTRO)
 376             {// Fake ALTRO
 377               fTriggerRawDigitMaker->Add( bunchlist );
 378             }//Fake ALTRO
 379         } // end while over channel
 380     } //end while over DDL's, of input stream
 381   fTriggerRawDigitMaker->PostProcess();
 382   TrimDigits(digitsArr);
 383 }
 384
 385
 386 void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr)
 387 {
 388   AliEMCALDigit *digit = 0;
 389   Int_t n = 0;
 390   Int_t nDigits = digitsArr->GetEntriesFast();
 391   TIter nextdigit(digitsArr);
 392   while ((digit = (AliEMCALDigit*) nextdigit())) {
 393     if (digit->GetType() == AliEMCALDigit::kLGnoHG) {
 394       AliDebug(1,Form("Remove digit with id %d, LGnoHG",digit->GetId()));
 395       digitsArr->Remove(digit);
 396     }
 397     else if(fTimeMin > digit->GetTime() || fTimeMax < digit->GetTime()) {
 398       digitsArr->Remove(digit);
 399       AliDebug(1,Form("Remove digit with id %d, Bad Time %e",digit->GetId(), digit->GetTime()));
 400     }
 401     else if (0 > digit->GetChi2()) {
 402       digitsArr->Remove(digit);
 403       AliDebug(1,Form("Remove digit with id %d, Bad Chi2 %e",digit->GetId(), digit->GetChi2()));
 404     }
 405     else {
 406       digit->SetIndexInList(n);
 407       n++;
 408     }
 409   }//while
 410
 411   digitsArr->Compress();
 412   AliDebug(1,Form("N Digits before trimming : %d; after array compression %d",nDigits,digitsArr->GetEntriesFast()));
 413 }
 414
 415
 416 Double_t
 417 AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
 418 {
 419   Double_t signal = 0.;
 420   Double_t tau    = par[2];
 421   Double_t n      = par[3];
 422   Double_t ped    = par[4];
 423   Double_t xx     = ( x[0] - par[1] + tau ) / tau ;
 424   if (xx <= 0)
 425     signal = ped ;
 426   else
 427     {
 428       signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
 429     }
 430   return signal ;
 431 }
 432
 433
 434 Bool_t AliEMCALRawUtils::RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH,
 435                                             Int_t * adcL, const Int_t keyErr) const
 436 {
 437   Bool_t lowGain = kFALSE ;
 438   TF1 signalF("signal", RawResponseFunction, 0, TIMEBINS, 5);
 439   signalF.SetParameter(0, damp) ;
 440   signalF.SetParameter(1, (dtime + fgTimeTrigger)/ TIMEBINWITH) ;
 441   signalF.SetParameter(2, TAU) ;
 442   signalF.SetParameter(3, ORDER);
 443   signalF.SetParameter(4, fgPedestalValue);
 444
 445   Double_t signal=0.0, noise=0.0;
 446   for (Int_t iTime = 0; iTime <  TIMEBINS; iTime++) {
 447     signal = signalF.Eval(iTime) ;
 448
 449     if(keyErr>0) {
 450       noise = gRandom->Gaus(0.,fgFEENoise);
 451       signal += noise;
 452     }
 453
 454     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
 455     if ( adcH[iTime] > MAXBINVALUE ){  // larger than 10 bits
 456       adcH[iTime] = MAXBINVALUE ;
 457       lowGain = kTRUE ;
 458     }
 459     signal /= HGLGFACTOR;
 460     adcL[iTime] =  static_cast<Int_t>(signal + 0.5) ;
 461     if ( adcL[iTime] > MAXBINVALUE )  // larger than 10 bits
 462       adcL[iTime] = MAXBINVALUE ;
 463   }
 464   return lowGain ;
 465 }
 466
 467
 468 //__________________________________________________________________
 469 Double_t AliEMCALRawUtils::RawResponseFunctionLog(Double_t *x, Double_t *par)
 470 {
 471   Double_t signal = 0. ;
 472   Double_t tau    = par[2];
 473   Double_t n      = par[3];
 474   Double_t xx     = ( x[0] - par[1] + tau ) / tau ;
 475
 476   if (xx < 0)
 477     {
 478       signal = par[0] - n*TMath::Log(TMath::Abs(xx)) + n * (1 - xx ) ;
 479     }
 480   else
 481     {
 482       signal = par[0] + n*TMath::Log(xx) + n * (1 - xx ) ;
 483     }
 484   return signal ;
 485 }
 486
 487
 488
 489 //__________________________________________________________________
 490 void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo)
 491 {
 492   // fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer(  Algo::kStandard );
 493   fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer( fitAlgo );
 494
 495   //fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer( kStandard );
 496
 497   fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done, for the new methods
 498   fRawAnalyzer->SetOverflowCut ( OVERFLOWCUT );
 499   fRawAnalyzer->SetAmpCut(fNoiseThreshold);
 500   fRawAnalyzer->SetFitArrayCut(fNoiseThreshold);
 501   //  return;
 502 }
 503
 504
 505