EMCAL/AliEMCALDigitizer.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 //
  20 //////////////////////////////////////////////////////////////////////////////
  21 // Class performs digitization of Summable digits from simulated data
  22 //
  23 // In addition it performs mixing/embedding of summable digits from different events.
  24 //
  25 // For each event 3 branches are created in TreeD:
  26 //   "EMCAL" - list of digits
  27 //   "EMCALTRG" - list of trigger digits
  28 //   "AliEMCALDigitizer" - AliEMCALDigitizer with all parameters used in digitization
  29 //
  30 //
  31 ////////////////////////////////////////////////////////////////////////////////////
  32 //
  33 //*-- Author: Sahal Yacoob (LBL)
  34 // based on : AliEMCALDigitizer
  35 // Modif:
  36 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
  37 //                           of new IO (a la PHOS)
  38 //  November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry
  39 //  July 2011 GCB: Digitizer modified to accomodate embedding.
  40 //                 Time calibration added. Decalibration possibility of energy and time added
  41 //_________________________________________________________________________________
  42
  43 // --- ROOT system ---
  44 #include <TROOT.h>
  45 #include <TTree.h>
  46 #include <TSystem.h>
  47 #include <TBenchmark.h>
  48 #include <TBrowser.h>
  49 #include <TObjectTable.h>
  50 #include <TRandom.h>
  51 #include <TF1.h>
  52 #include <cassert>
  53
  54 // --- AliRoot header files ---
  55 #include "AliLog.h"
  56 #include "AliRun.h"
  57 #include "AliDigitizationInput.h"
  58 #include "AliRunLoader.h"
  59 #include "AliCDBManager.h"
  60 #include "AliCDBEntry.h"
  61 #include "AliEMCALDigit.h"
  62 #include "AliEMCAL.h"
  63 #include "AliEMCALLoader.h"
  64 #include "AliEMCALDigitizer.h"
  65 #include "AliEMCALSDigitizer.h"
  66 #include "AliEMCALGeometry.h"
  67 #include "AliEMCALCalibData.h"
  68 #include "AliEMCALSimParam.h"
  69 #include "AliEMCALRawDigit.h"
  70 #include "AliCaloCalibPedestal.h"
  71
  72   namespace
  73   {
  74     Double_t HeavisideTheta(Double_t x)
  75     {
  76       Double_t signal = 0.;
  77
  78       if (x > 0.) signal = 1.;
  79
  80       return signal;
  81     }
  82
  83     Double_t AnalogFastORFunction(Double_t *x, Double_t *par)
  84     {
  85       Double_t v0 = par[0];
  86       Double_t t0 = par[1];
  87       Double_t tr = par[2];
  88
  89       Double_t R1 = 1000.;
  90       Double_t C1 = 33e-12;
  91       Double_t R2 = 1800;
  92       Double_t C2 = 22e-12;
  93
  94       Double_t t  =   x[0];
  95
  96       return (((0.8*(-((TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-t + t0)/(C1*R1))*
  97                         TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2)) +
  98                      C1*C2*R1*R2*(1 - (C2*TMath::Power(TMath::E(),(-t + t0)/(C2*R2))*R2)/(-(C1*R1) + C2*R2)))*v0*
  99                 HeavisideTheta(t - t0))/tr
 100                - (0.8*(C1*C2*R1*R2 -
 101                        (TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C1*R1))*
 102                         TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2) +
 103                        (C1*TMath::Power(C2,2)*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C2*R2))*
 104                         R1*TMath::Power(R2,2))/(C1*R1 - C2*R2))*v0*
 105                   HeavisideTheta(t - t0 - 1.25*tr))/tr)/(C2*R1));
 106     }
 107   }
 108
 109 ClassImp(AliEMCALDigitizer)
 110
 111
 112 //____________________________________________________________________________
 113 AliEMCALDigitizer::AliEMCALDigitizer()
 114   : AliDigitizer("",""),
 115     fDefaultInit(kTRUE),
 116     fDigitsInRun(0),
 117     fInit(0),
 118     fInput(0),
 119     fInputFileNames(0x0),
 120     fEventNames(0x0),
 121     fDigitThreshold(0),
 122     fMeanPhotonElectron(0),
 123     fGainFluctuations(0),
 124     fPinNoise(0),
 125     fTimeNoise(0),
 126     fTimeDelay(0),
 127     fTimeResolutionPar0(0),
 128     fTimeResolutionPar1(0),
 129     fADCchannelEC(0),
 130     fADCpedestalEC(0),
 131     fADCchannelECDecal(0),
 132     fTimeChannel(0),
 133     fTimeChannelDecal(0),
 134     fNADCEC(0),
 135     fEventFolderName(""),
 136     fFirstEvent(0),
 137     fLastEvent(0),
 138     fCalibData(0x0),
 139     fSDigitizer(0x0)
 140 {
 141   // ctor
 142   InitParameters() ;
 143   fDigInput = 0 ;                     // We work in the standalone mode
 144 }
 145
 146 //____________________________________________________________________________
 147 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
 148   : AliDigitizer("EMCALDigitizer", alirunFileName),
 149     fDefaultInit(kFALSE),
 150     fDigitsInRun(0),
 151     fInit(0),
 152     fInput(0),
 153     fInputFileNames(0),
 154     fEventNames(0),
 155     fDigitThreshold(0),
 156     fMeanPhotonElectron(0),
 157     fGainFluctuations(0),
 158     fPinNoise(0),
 159     fTimeNoise(0),
 160     fTimeDelay(0),
 161     fTimeResolutionPar0(0),
 162     fTimeResolutionPar1(0),
 163     fADCchannelEC(0),
 164     fADCpedestalEC(0),
 165     fADCchannelECDecal(0),
 166     fTimeChannel(0),
 167     fTimeChannelDecal(0),
 168     fNADCEC(0),
 169     fEventFolderName(eventFolderName),
 170     fFirstEvent(0),
 171     fLastEvent(0),
 172     fCalibData(0x0),
 173     fSDigitizer(0x0)
 174 {
 175   // ctor
 176   InitParameters() ;
 177   Init() ;
 178   fDigInput = 0 ;                     // We work in the standalone mode
 179 }
 180
 181 //____________________________________________________________________________
 182 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
 183   : AliDigitizer(d.GetName(),d.GetTitle()),
 184     fDefaultInit(d.fDefaultInit),
 185     fDigitsInRun(d.fDigitsInRun),
 186     fInit(d.fInit),
 187     fInput(d.fInput),
 188     fInputFileNames(d.fInputFileNames),
 189     fEventNames(d.fEventNames),
 190     fDigitThreshold(d.fDigitThreshold),
 191     fMeanPhotonElectron(d.fMeanPhotonElectron),
 192     fGainFluctuations(d.fGainFluctuations),
 193     fPinNoise(d.fPinNoise),
 194     fTimeNoise(d.fTimeNoise),
 195     fTimeDelay(d.fTimeDelay),
 196     fTimeResolutionPar0(d.fTimeResolutionPar0),
 197     fTimeResolutionPar1(d.fTimeResolutionPar1),
 198     fADCchannelEC(d.fADCchannelEC),
 199     fADCpedestalEC(d.fADCpedestalEC),
 200     fADCchannelECDecal(d.fADCchannelECDecal),
 201     fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
 202     fNADCEC(d.fNADCEC),
 203     fEventFolderName(d.fEventFolderName),
 204     fFirstEvent(d.fFirstEvent),
 205     fLastEvent(d.fLastEvent),
 206     fCalibData(d.fCalibData),
 207     fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
 208 {
 209   // copyy ctor
 210 }
 211
 212 //____________________________________________________________________________
 213 AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
 214   : AliDigitizer(rd,"EMCALDigitizer"),
 215     fDefaultInit(kFALSE),
 216     fDigitsInRun(0),
 217     fInit(0),
 218     fInput(0),
 219     fInputFileNames(0),
 220     fEventNames(0),
 221     fDigitThreshold(0),
 222     fMeanPhotonElectron(0),
 223     fGainFluctuations(0),
 224     fPinNoise(0.),
 225     fTimeNoise(0.),
 226     fTimeDelay(0.),
 227     fTimeResolutionPar0(0.),
 228     fTimeResolutionPar1(0.),
 229     fADCchannelEC(0),
 230     fADCpedestalEC(0),
 231     fADCchannelECDecal(0),
 232     fTimeChannel(0),
 233     fTimeChannelDecal(0),
 234     fNADCEC(0),
 235     fEventFolderName(0),
 236     fFirstEvent(0),
 237     fLastEvent(0),
 238     fCalibData(0x0),
 239     fSDigitizer(0x0)
 240 {
 241   // ctor Init() is called by RunDigitizer
 242   fDigInput = rd ;
 243   fEventFolderName = fDigInput->GetInputFolderName(0) ;
 244   SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
 245   InitParameters() ;
 246 }
 247
 248 //____________________________________________________________________________
 249   AliEMCALDigitizer::~AliEMCALDigitizer()
 250 {
 251   //dtor
 252   delete [] fInputFileNames ;
 253   delete [] fEventNames ;
 254   if (fSDigitizer) delete fSDigitizer;
 255 }
 256
 257 //____________________________________________________________________________
 258 void AliEMCALDigitizer::Digitize(Int_t event)
 259 {
 260   // Makes the digitization of the collected summable digits
 261   // for this it first creates the array of all EMCAL modules
 262   // filled with noise and after that adds contributions from
 263   // SDigits. This design helps to avoid scanning over the
 264   // list of digits to add  contribution of any new SDigit.
 265   //
 266   // JLK 26-Jun-2008
 267   // Note that SDigit energy info is stored as an amplitude, so we
 268   // must call the Calibrate() method of the SDigitizer to convert it
 269   // back to an energy in GeV before adding it to the Digit
 270   //
 271
 272   AliRunLoader *rl = AliRunLoader::Instance();
 273   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 274
 275   if(!emcalLoader)
 276   {
 277     AliFatal("EMCALLoader is NULL!") ;
 278     return; // not needed but in case coverity complains ...
 279   }
 280
 281   Int_t readEvent = event ;
 282   if (fDigInput) // fDigInput is data member from AliDigitizer
 283     readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
 284   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 285                   readEvent, GetTitle(), fEventFolderName.Data())) ;
 286   rl->GetEvent(readEvent);
 287
 288   TClonesArray * digits = emcalLoader->Digits() ;
 289   digits->Delete() ;  //correct way to clear array when memory is
 290   //allocated by objects stored in array
 291
 292   // Load Geometry
 293   if (!rl->GetAliRun())
 294   {
 295     AliFatal("Could not get AliRun from runLoader");
 296     return; // not needed but in case coverity complains ...
 297   }
 298
 299   AliEMCAL * emcal  = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
 300   AliEMCALGeometry *geom = emcal->GetGeometry();
 301   static int nEMC = geom->GetNCells();//max number of digits possible
 302   AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
 303
 304   digits->Expand(nEMC) ;
 305
 306   // RS create a digitizer on the fly
 307   if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
 308   fSDigitizer->SetEventRange(0, -1) ;
 309
 310   //-------------------------------------------------------
 311   //take all the inputs to add together and load the SDigits
 312   TObjArray * sdigArray = new TObjArray(fInput) ;
 313   sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
 314
 315   Int_t i ;
 316   Int_t embed = kFALSE;
 317   for(i = 1 ; i < fInput ; i++)
 318   {
 319     TString tempo(fEventNames[i]) ;
 320     tempo += i ;
 321
 322     AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ;
 323
 324     if (!rl2)
 325       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
 326
 327     if(!rl2)
 328     {
 329       AliFatal("Run Loader of second event not available!");
 330       return; // not needed but in case coverity complains ...
 331     }
 332
 333     if (fDigInput)
 334       readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
 335
 336     Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
 337
 338     rl2->LoadSDigits();
 339     //     rl2->LoadDigits();
 340     rl2->GetEvent(readEvent);
 341
 342     if(!rl2->GetDetectorLoader("EMCAL"))
 343     {
 344       AliFatal("Loader of second input not found");
 345       return; // not needed but in case coverity complains ...
 346     }
 347
 348     AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
 349
 350     if(!emcalLoader2)
 351     {
 352       AliFatal("EMCAL Loader of second event not available!");
 353       return; // not needed but in case coverity complains ...
 354     }
 355
 356     //Merge 2 simulated sdigits
 357     if(!emcalLoader2->SDigits()) continue;
 358
 359     TClonesArray* sdigits2 = emcalLoader2->SDigits();
 360     sdigArray->AddAt(sdigits2, i) ;
 361
 362     // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
 363     // do not smear energy of embedded, do not add noise to any sdigits
 364     if( sdigits2->GetEntriesFast() <= 0 ) continue;
 365
 366     //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
 367     AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
 368     if( digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded ) embed = kTRUE;
 369
 370   }// input loop
 371
 372   //--------------------------------
 373   //Find the first tower with signal
 374   Int_t nextSig = nEMC + 1 ;
 375   TClonesArray * sdigits ;
 376   for(i = 0 ; i < fInput ; i++)
 377   {
 378     if(i > 0 && embed) continue;
 379
 380     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 381     if(!sdigits)
 382     {
 383       AliDebug(1,"Null sdigit pointer");
 384       continue;
 385     }
 386
 387     if (sdigits->GetEntriesFast() <= 0 )
 388     {
 389       AliDebug(1, "No sdigits entries");
 390       continue;
 391     }
 392
 393     AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
 394     if(!sd)
 395     {
 396       AliDebug(1, "NULL sdigit pointer");
 397       continue;
 398     }
 399
 400     Int_t curNext = sd->GetId() ;
 401     if(curNext < nextSig)
 402     nextSig = curNext ;
 403     AliDebug(1,Form("input %i : #sdigits %i \n",i, sdigits->GetEntriesFast()));
 404
 405   }// input loop
 406
 407   AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
 408
 409   TArrayI index(fInput) ;
 410   index.Reset() ;  //Set all indexes to zero
 411
 412   AliEMCALDigit * digit ;
 413   AliEMCALDigit * curSDigit ;
 414
 415   //---------------------------------------------
 416   //Put Noise contribution, smear time and energy
 417   Float_t timeResolution = 0;
 418   Int_t absID = -1 ;
 419   for(absID = 0; absID < nEMC; absID++)
 420   { // Nov 30, 2006 by PAI; was from 1 to nEMC
 421
 422     if (IsDead(absID)) continue; // Don't instantiate dead digits
 423
 424     Float_t energy = 0 ;
 425
 426     // amplitude set to zero, noise will be added later
 427     Float_t noiseTime = 0.;
 428     if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
 429     new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
 430     //look if we have to add signal?
 431     digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
 432
 433     if (!digit)
 434     {
 435       AliDebug(1,"Digit pointer is null");
 436       continue;
 437     }
 438
 439     if(absID==nextSig)
 440     {
 441       // Calculate time as time of the largest digit
 442       Float_t time = digit->GetTime() ;
 443       Float_t aTime= digit->GetAmplitude() ;
 444
 445       // loop over input
 446       Int_t nInputs = fInput;
 447       if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
 448       for(i = 0; i< nInputs ; i++)
 449       {
 450         //loop over (possible) merge sources
 451         TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
 452         if(sdtclarr)
 453         {
 454           Int_t sDigitEntries = sdtclarr->GetEntriesFast();
 455
 456           if(sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
 457           else                          curSDigit = 0 ;
 458
 459           //May be several digits will contribute from the same input
 460           while(curSDigit && (curSDigit->GetId() == absID))
 461           {
 462             //Shift primary to separate primaries belonging different inputs
 463             Int_t primaryoffset = i ;
 464             if(fDigInput) primaryoffset = fDigInput->GetMask(i) ;
 465             curSDigit->ShiftPrimary(primaryoffset) ;
 466
 467             if(curSDigit->GetAmplitude()>aTime)
 468             {
 469               aTime = curSDigit->GetAmplitude();
 470               time  = curSDigit->GetTime();
 471             }
 472
 473             *digit = *digit + *curSDigit ;  //adds amplitudes of each digit
 474
 475             index[i]++ ;
 476
 477             if( sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
 478             else                           curSDigit = 0 ;
 479           }// while
 480         }// source exists
 481       }// loop over merging sources
 482
 483       //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
 484       energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
 485
 486       // add fluctuations for photo-electron creation
 487       // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
 488       Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
 489       energy       *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
 490
 491       //calculate and set time
 492       digit->SetTime(time) ;
 493
 494       //Find next signal module
 495       nextSig = nEMC + 1 ;
 496       for(i = 0 ; i < nInputs ; i++)
 497       {
 498         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 499
 500         if(sdigits)
 501         {
 502           Int_t curNext = nextSig ;
 503           if(sdigits->GetEntriesFast() > index[i])
 504           {
 505             AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
 506             if ( tmpdigit ) curNext = tmpdigit->GetId() ;
 507           }
 508
 509           if(curNext < nextSig) nextSig = curNext ;
 510         }// sdigits exist
 511       } // input loop
 512
 513     }//absID==nextSig
 514
 515     // add the noise now, no need for embedded events with real data
 516     if(!embed)
 517       energy += gRandom->Gaus(0., fPinNoise) ;
 518
 519     // JLK 26-June-2008
 520     //Now digitize the energy according to the fSDigitizer method,
 521     //which merely converts the energy to an integer.  Later we will
 522     //check that the stored value matches our allowed dynamic ranges
 523     digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
 524
 525     //Set time resolution with final energy
 526     timeResolution = GetTimeResolution(energy);
 527     digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
 528     AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
 529                      absID, energy, nextSig));
 530     //Add delay to time
 531     digit->SetTime(digit->GetTime()+fTimeDelay) ;
 532
 533   } // for(absID = 0; absID < nEMC; absID++)
 534
 535   //---------------------------------------------------------
 536   //Embed simulated amplitude (and time?) to real data digits
 537   if ( embed )
 538   {
 539     for(Int_t input = 1; input<fInput; input++)
 540     {
 541       TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
 542       if(!realDigits)
 543       {
 544         AliDebug(1,"No sdigits to merge\n");
 545         continue;
 546       }
 547
 548       for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
 549       {
 550         AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
 551         if ( !digit2 ) continue;
 552
 553         digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
 554         if ( !digit ) continue;
 555
 556         // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
 557         // and multiply to get a big int.
 558         Float_t time2 = digit2->GetTime();
 559         Float_t e2    = digit2->GetAmplitude();
 560         CalibrateADCTime(e2,time2,digit2->GetId());
 561         Float_t amp2  = fSDigitizer->Digitize(e2);
 562
 563         Float_t e0    = digit ->GetAmplitude();
 564         if(e0>0.01)
 565           AliDebug(1,Form("digit 1: Abs ID %d, amp %f, type %d, time %e; digit2: Abs ID %d, amp %f, type %d, time %e\n",
 566                           digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
 567                           digit2->GetId(),amp2,                   digit2->GetType(), time2           ));
 568
 569         // Sum amplitudes, change digit amplitude
 570         digit->SetAmplitude( digit->GetAmplitude() + amp2);
 571
 572         // Rough assumption, assign time of the largest of the 2 digits,
 573         // data or signal to the final digit.
 574         if(amp2 > digit->GetAmplitude())  digit->SetTime(time2);
 575
 576         //Mark the new digit as embedded
 577         digit->SetType(AliEMCALDigit::kEmbedded);
 578
 579         if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
 580           AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
 581                           digit->GetId(), digit->GetAmplitude(), digit->GetType()));
 582       }//loop digit 2
 583     }//input loop
 584   }//embed
 585
 586   //JLK is it better to call Clear() here?
 587   delete sdigArray ; //We should not delete its contents
 588
 589   //------------------------------
 590   //remove digits below thresholds
 591   // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
 592   // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
 593   // before merge in the same loop real data digits if available
 594   Float_t energy = 0;
 595   Float_t time   = 0;
 596   for(i = 0 ; i < nEMC ; i++)
 597   {
 598     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
 599     if ( !digit ) continue;
 600
 601     //Then get the energy in GeV units.
 602     energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
 603
 604     //Then digitize using the calibration constants of the ocdb
 605     Float_t ampADC = energy;
 606     DigitizeEnergyTime(ampADC, time, digit->GetId())  ;
 607
 608     if(ampADC < fDigitThreshold)
 609       digits->RemoveAt(i) ;
 610
 611   } // digit loop
 612
 613   digits->Compress() ;
 614
 615   Int_t ndigits = digits->GetEntriesFast() ;
 616
 617   //---------------------------------------------------------------
 618   //JLK 26-June-2008
 619   //After we have done the summing and digitizing to create the
 620   //digits, now we want to calibrate the resulting amplitude to match
 621   //the dynamic range of our real data.
 622   for (i = 0 ; i < ndigits ; i++)
 623   {
 624     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
 625     if( !digit ) continue ;
 626
 627     digit->SetIndexInList(i) ;
 628
 629     time   = digit->GetTime();
 630     digit->SetTime(time);
 631
 632     energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
 633
 634     Float_t ampADC = energy;
 635     DigitizeEnergyTime(ampADC, time, digit->GetId());
 636
 637     digit->SetAmplitude(ampADC) ;
 638     // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
 639   }//Digit loop
 640
 641 }
 642
 643 //_____________________________________________________________________
 644 void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
 645 {
 646   // JLK 26-June-2008
 647   // Returns digitized value of the energy in a cell absId
 648   // using the calibration constants stored in the OCDB
 649   // or default values if no CalibData object is found.
 650   // This effectively converts everything to match the dynamic range
 651   // of the real data we will collect
 652   //
 653   // Load Geometry
 654   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 655
 656   if (geom==0)
 657   {
 658     AliFatal("Did not get geometry from EMCALLoader");
 659     return;
 660   }
 661
 662   Int_t iSupMod = -1;
 663   Int_t nModule = -1;
 664   Int_t nIphi   = -1;
 665   Int_t nIeta   = -1;
 666   Int_t iphi    = -1;
 667   Int_t ieta    = -1;
 668   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 669   if(!bCell)
 670   Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
 671   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 672
 673   if(fCalibData)
 674   {
 675     fADCpedestalEC     = fCalibData->GetADCpedestal     (iSupMod,ieta,iphi);
 676     fADCchannelEC      = fCalibData->GetADCchannel      (iSupMod,ieta,iphi);
 677     fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
 678     fTimeChannel       = fCalibData->GetTimeChannel     (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
 679     fTimeChannelDecal  = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
 680   }
 681
 682   //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
 683   energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal   ;
 684   time  += fTimeChannel-fTimeChannelDecal;
 685
 686   if ( energy > fNADCEC ) energy =  fNADCEC ;
 687 }
 688
 689 //_____________________________________________________________________
 690 void AliEMCALDigitizer::DecalibrateTrigger(AliEMCALDigit *digit)
 691 {
 692   // Decalibrate, used in Trigger digits
 693
 694   if ( !fCalibData ) return ;
 695
 696   // Load Geometry
 697   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 698
 699   if (!geom)
 700   {
 701     AliFatal("Did not get geometry from EMCALLoader");
 702     return;
 703   }
 704
 705   // Recover the digit location in row-column-SM
 706   Int_t absId   = digit->GetId();
 707   Int_t iSupMod = -1;
 708   Int_t nModule = -1;
 709   Int_t nIphi   = -1;
 710   Int_t nIeta   = -1;
 711   Int_t iphi    = -1;
 712   Int_t ieta    = -1;
 713
 714   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 715   if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
 716
 717   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 718
 719   // Now decalibrate
 720   Float_t adcChannel       = fCalibData->GetADCchannel      (iSupMod,ieta,iphi);
 721   Float_t adcChannelOnline = fCalibData->GetADCchannelOnline(iSupMod,ieta,iphi);
 722   //printf("calib param for (SM%d,col%d,row%d): %1.4f - online param: %1.4f \n",iSupMod,ieta,iphi,adcChannel,adcChannelOnline);
 723   Float_t factor = 1./(adcChannel/adcChannelOnline);
 724
 725   *digit = *digit * factor;
 726
 727 }
 728
 729 //_____________________________________________________________________
 730 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
 731 {
 732   // Returns the energy in a cell absId with a given adc value
 733   // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
 734   // Used in case of embedding, transform ADC counts from real event into energy
 735   // so that we can add the energy of the simulated sdigits which are in energy
 736   // units.
 737   // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
 738   // or with time out of window
 739
 740   // Load Geometry
 741   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 742
 743   if (!geom)
 744   {
 745     AliFatal("Did not get geometry from EMCALLoader");
 746     return;
 747   }
 748
 749   Int_t iSupMod = -1;
 750   Int_t nModule = -1;
 751   Int_t nIphi   = -1;
 752   Int_t nIeta   = -1;
 753   Int_t iphi    = -1;
 754   Int_t ieta    = -1;
 755   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 756   if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
 757   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 758
 759   if(fCalibData)
 760   {
 761     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 762     fADCchannelEC  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 763     fTimeChannel   = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
 764   }
 765
 766   adc   = adc * fADCchannelEC - fADCpedestalEC;
 767   time -= fTimeChannel;
 768
 769 }
 770
 771
 772 //____________________________________________________________________________
 773 void AliEMCALDigitizer::Digitize(Option_t *option)
 774 {
 775   // Steering method to process digitization for events
 776   // in the range from fFirstEvent to fLastEvent.
 777   // This range is optionally set by SetEventRange().
 778   // if fLastEvent=-1, then process events until the end.
 779   // by default fLastEvent = fFirstEvent (process only one event)
 780
 781   if (!fInit)
 782   { // to prevent overwrite existing file
 783     Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
 784     return ;
 785   }
 786
 787   if (strstr(option,"print"))
 788   {
 789     Print();
 790     return ;
 791   }
 792
 793   if(strstr(option,"tim"))
 794     gBenchmark->Start("EMCALDigitizer");
 795
 796   AliRunLoader *rl = AliRunLoader::Instance();
 797
 798   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 799   if(!emcalLoader)
 800   {
 801     AliFatal("Did not get the  Loader");
 802     return; // coverity
 803   }
 804
 805   if (fLastEvent == -1)
 806     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 807   else if (fDigInput)
 808     fLastEvent = fFirstEvent ; // what is this ??
 809
 810   Int_t nEvents = fLastEvent - fFirstEvent + 1;
 811   Int_t ievent  = -1;
 812
 813   AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
 814   if(!emcal)
 815   {
 816     AliFatal("Did not get the AliEMCAL pointer");
 817     return; // coverity
 818   }
 819
 820   AliEMCALGeometry *geom = emcal->GetGeometry();
 821   if(!geom)
 822   {
 823     AliFatal("Geometry pointer null");
 824     return; // fix for coverity
 825   }
 826
 827   const Int_t nTRU = geom->GetNTotalTRU();
 828   TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit",    nTRU*96);
 829   TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
 830
 831   rl->LoadSDigits("EMCAL");
 832   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++)
 833   {
 834     rl->GetEvent(ievent);
 835
 836     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 837
 838     WriteDigits() ;
 839
 840     //Trigger Digits
 841     //-------------------------------------
 842
 843     Digits2FastOR(digitsTMP, digitsTRG);
 844
 845     WriteDigits(digitsTRG);
 846
 847     (emcalLoader->TreeD())->Fill();
 848
 849     emcalLoader->WriteDigits(   "OVERWRITE");
 850
 851     Unload();
 852
 853     digitsTRG  ->Delete();
 854     digitsTMP  ->Delete();
 855
 856     //-------------------------------------
 857
 858     if(strstr(option,"deb"))
 859       PrintDigits(option);
 860     if(strstr(option,"table")) gObjectTable->Print();
 861
 862     //increment the total number of Digits per run
 863     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
 864   }//loop
 865
 866   if(strstr(option,"tim"))
 867   {
 868     gBenchmark->Stop("EMCALDigitizer");
 869     Float_t cputime   = gBenchmark->GetCpuTime("EMCALDigitizer");
 870     Float_t avcputime = cputime;
 871     if(nEvents==0) avcputime = 0 ;
 872     AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
 873   }
 874 }
 875
 876 //__________________________________________________________________
 877 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
 878 {
 879   // Assign a smeared time to the digit, from observed distribution in LED system (?).
 880   // From F. Blanco
 881   Float_t res = -1;
 882   if (energy > 0)
 883   {
 884     res = TMath::Sqrt(fTimeResolutionPar0 +
 885                       fTimeResolutionPar1/(energy*energy) );
 886   }
 887   // parametrization above is for ns. Convert to seconds:
 888   return res*1e-9;
 889 }
 890
 891 //____________________________________________________________________________
 892 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
 893 {
 894   // FEE digits afterburner to produce TRG digits
 895   // we are only interested in the FEE digit deposited energy
 896   // to be converted later into a voltage value
 897
 898   // push the FEE digit to its associated FastOR (numbered from 0:95)
 899   // TRU is in charge of summing module digits
 900
 901   AliRunLoader *runLoader = AliRunLoader::Instance();
 902
 903   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
 904   if (!emcalLoader) AliFatal("Did not get the  Loader");
 905
 906   const  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 907   if(!geom)
 908   {
 909     AliFatal("Geometry pointer null");
 910     return; // fix for coverity
 911   }
 912
 913   // build FOR from simulated digits
 914   // and xfer to the corresponding TRU input (mapping)
 915
 916   TClonesArray* sdigits = emcalLoader->SDigits();
 917
 918   TClonesArray *digits = (TClonesArray*)sdigits->Clone();
 919
 920   AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
 921
 922   TIter NextDigit(digits);
 923   while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
 924   {
 925     if (IsDead(digit)) continue;
 926
 927     DecalibrateTrigger(digit);
 928
 929     Int_t id = digit->GetId();
 930
 931     Int_t trgid;
 932     if (geom->GetFastORIndexFromCellIndex(id , trgid))
 933     {
 934       AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
 935
 936       AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
 937
 938       if (!d)
 939       {
 940         new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
 941         d = (AliEMCALDigit*)digitsTMP->At(trgid);
 942         d->SetId(trgid);
 943       }
 944       else
 945       {
 946         *d = *d + *digit;
 947       }
 948     }
 949   }
 950
 951   if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
 952
 953   Int_t     nSamples = geom->GetNTotalTRU();
 954   Int_t *timeSamples = new Int_t[nSamples];
 955
 956   NextDigit = TIter(digitsTMP);
 957   while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
 958   {
 959     if (digit)
 960     {
 961       Int_t     id = digit->GetId();
 962       Float_t time = 50.e-9;
 963
 964       Double_t depositedEnergy = 0.;
 965       for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
 966
 967       if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
 968
 969       // FIXME: Check digit time!
 970       if (depositedEnergy) {
 971         depositedEnergy += gRandom->Gaus(0., .08);
 972         DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
 973
 974         for (Int_t j=0;j<nSamples;j++) {
 975           if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
 976           timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
 977         }
 978
 979         new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
 980
 981         if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
 982       }
 983     }
 984   }
 985
 986   delete [] timeSamples;
 987
 988   if (digits && digits->GetEntriesFast()) digits->Delete();
 989 }
 990
 991 //____________________________________________________________________________
 992 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
 993 {
 994   // parameters:
 995   // id: 0..95
 996   const Int_t    reso = 12;      // 11-bit resolution ADC
 997   const Double_t vFSR = 2.;      // Full scale input voltage range 2V (p-p)
 998 //const Double_t dNe  = 125;     // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
 999   const Double_t dNe  = 125/1.3; // F-ALTRO max V. FEE: factor 4
1000   const Double_t vA   = .136e-6; // CSP output range: 0.136uV/e-
1001   const Double_t rise = 50e-9;   // rise time (10-90%) of the FastOR signal before shaping
1002
1003   const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
1004
1005   Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
1006
1007   TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
1008   signalF.SetParameter( 0,   vV );
1009   signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
1010   signalF.SetParameter( 2, rise );
1011
1012   for (Int_t iTime=0; iTime<nSamples; iTime++)
1013   {
1014     // FIXME: add noise (probably not simply Gaussian) according to DA measurements
1015     // probably plan an access to OCDB
1016     Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
1017     if (TMath::Abs(sig) > vFSR/2.) {
1018       AliError("Signal overflow!");
1019       timeSamples[iTime] = (1 << reso) - 1;
1020     } else {
1021       AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
1022       timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
1023     }
1024   }
1025 }
1026
1027 //____________________________________________________________________________
1028 Bool_t AliEMCALDigitizer::Init()
1029 {
1030   // Makes all memory allocations
1031   fInit = kTRUE ;
1032   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1033
1034   if ( emcalLoader == 0 ) {
1035     Fatal("Init", "Could not obtain the AliEMCALLoader");
1036     return kFALSE;
1037   }
1038
1039   fFirstEvent = 0 ;
1040   fLastEvent = fFirstEvent ;
1041
1042   if(fDigInput)
1043     fInput = fDigInput->GetNinputs() ;
1044   else
1045     fInput           = 1 ;
1046
1047   fInputFileNames  = new TString[fInput] ;
1048   fEventNames      = new TString[fInput] ;
1049   fInputFileNames[0] = GetTitle() ;
1050   fEventNames[0]     = fEventFolderName.Data() ;
1051   Int_t index ;
1052   for (index = 1 ; index < fInput ; index++) {
1053     fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1054     TString tempo = fDigInput->GetInputFolderName(index) ;
1055     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
1056   }
1057
1058   //Calibration instance
1059   fCalibData = emcalLoader->CalibData();
1060   return fInit ;
1061 }
1062
1063 //____________________________________________________________________________
1064 void AliEMCALDigitizer::InitParameters()
1065 {
1066   // Parameter initialization for digitizer
1067
1068   // Get the parameters from the OCDB via the loader
1069   AliRunLoader *rl = AliRunLoader::Instance();
1070   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1071   AliEMCALSimParam * simParam = 0x0;
1072   if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1073
1074   if(!simParam){
1075     simParam = AliEMCALSimParam::GetInstance();
1076     AliWarning("Simulation Parameters not available in OCDB?");
1077   }
1078
1079   fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400;  // electrons per GeV
1080   fGainFluctuations   = simParam->GetGainFluctuations()   ; // 15.0;
1081
1082   fPinNoise           = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1083   if (fPinNoise < 0.0001 )
1084     Warning("InitParameters", "No noise added\n") ;
1085   fTimeNoise          = simParam->GetTimeNoise(); // 1.28E-5 s
1086   fDigitThreshold     = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1087   fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1088   fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1089   fTimeDelay          = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1090
1091   // These defaults are normally not used.
1092   // Values are read from calibration database instead
1093   fADCchannelEC       = 0.0162; // Nominal value set online for most of the channels, MIP peak sitting at ~266./16/1024
1094   fADCpedestalEC      = 0.0 ;   // GeV
1095   fADCchannelECDecal  = 1.0;    // No decalibration by default
1096   fTimeChannel        = 0.0;    // No time calibration by default
1097   fTimeChannelDecal   = 0.0;    // No time decalibration by default
1098
1099   fNADCEC             = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
1100
1101   AliDebug(2,Form("Mean Photon Electron %d, Gain Fluct. %2.1f; Noise: APD %f, Time %f; Digit Threshold %d,Time Resolution Par0 %g Par1 %g,NADCEC %d",
1102                   fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1103
1104 }
1105
1106 //__________________________________________________________________
1107 void AliEMCALDigitizer::Print1(Option_t * option)
1108 { // 19-nov-04 - just for convenience
1109   Print();
1110   PrintDigits(option);
1111 }
1112
1113 //__________________________________________________________________
1114 void AliEMCALDigitizer::Print (Option_t * ) const
1115 {
1116   // Print Digitizer's parameters
1117   printf("Print: \n------------------- %s -------------", GetName() ) ;
1118   if( strcmp(fEventFolderName.Data(), "") != 0 ){
1119     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
1120
1121     Int_t nStreams ;
1122     if (fDigInput)
1123       nStreams =  GetNInputStreams() ;
1124     else
1125       nStreams = fInput ;
1126
1127     AliRunLoader *rl=0;
1128
1129     Int_t index = 0 ;
1130     for (index = 0 ; index < nStreams ; index++) {
1131       TString tempo(fEventNames[index]) ;
1132       tempo += index ;
1133       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1134         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1135       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1136       if(emcalLoader){
1137         TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1138         if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1139           fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
1140         printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1141       }//loader
1142     }
1143
1144     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1145
1146     if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1147     else printf("\nNULL LOADER");
1148
1149     printf("\nWith following parameters:\n") ;
1150     printf("    Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1151     printf("    Threshold  in Tower  (fDigitThreshold) = %d\n", fDigitThreshold)  ;
1152     printf("---------------------------------------------------\n")  ;
1153   }
1154   else
1155     printf("Print: AliEMCALDigitizer not initialized") ;
1156 }
1157
1158 //__________________________________________________________________
1159 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1160 {
1161   //utility method for printing digit information
1162
1163   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1164   if(emcalLoader){
1165     TClonesArray * digits  = emcalLoader->Digits() ;
1166     TClonesArray * sdigits = emcalLoader->SDigits() ;
1167
1168     printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1169     printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1170
1171     if(strstr(option,"all")){
1172       //loop over digits
1173       AliEMCALDigit * digit;
1174       printf("\nEMC digits (with primaries):\n")  ;
1175       printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
1176       Int_t index ;
1177       for (index = 0 ; index < digits->GetEntries() ; index++) {
1178         digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1179         if(digit){
1180           printf("\n%6d  %8f    %6.5e %4d      %2d : ",
1181                  digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1182           Int_t iprimary;
1183           for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1184             printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1185           }
1186         }// digit exists
1187       }// loop
1188     }
1189     printf("\n");
1190   }// loader exists
1191   else printf("NULL LOADER, cannot print\n");
1192 }
1193
1194 //__________________________________________________________________
1195 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1196 {
1197   // Calculates the time signal generated by noise
1198   //printf("Time noise %e\n",fTimeNoise);
1199   return gRandom->Rndm() * fTimeNoise;
1200 }
1201
1202 //__________________________________________________________________
1203 void AliEMCALDigitizer::Unload()
1204 {
1205   // Unloads the SDigits and Digits
1206   AliRunLoader *rl=0;
1207
1208   Int_t i ;
1209   for(i = 1 ; i < fInput ; i++){
1210     TString tempo(fEventNames[i]) ;
1211     tempo += i;
1212     if ((rl = AliRunLoader::GetRunLoader(tempo)))
1213       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1214   }
1215   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1216   if(emcalLoader)emcalLoader->UnloadDigits() ;
1217   else AliFatal("Did not get the loader");
1218 }
1219
1220 //_________________________________________________________________________________________
1221 void AliEMCALDigitizer::WriteDigits()
1222 { // Makes TreeD in the output file.
1223   // Check if branch already exists:
1224   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1225   //      already existing branches.
1226   //   else creates branch with Digits, named "EMCAL", title "...",
1227   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1228   //      and names of files, from which digits are made.
1229
1230   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1231
1232   if(emcalLoader){
1233     const TClonesArray * digits = emcalLoader->Digits() ;
1234     TTree * treeD = emcalLoader->TreeD();
1235     if ( !treeD ) {
1236       emcalLoader->MakeDigitsContainer();
1237       treeD = emcalLoader->TreeD();
1238     }
1239
1240     // -- create Digits branch
1241     Int_t bufferSize = 32000 ;
1242     TBranch * digitsBranch = 0;
1243     if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1244       digitsBranch->SetAddress(&digits);
1245       AliWarning("Digits Branch already exists. Not all digits will be visible");
1246     }
1247     else
1248       treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1249     //digitsBranch->SetTitle(fEventFolderName);
1250
1251     //  treeD->Fill() ;
1252     /*
1253      emcalLoader->WriteDigits("OVERWRITE");
1254      emcalLoader->WriteDigitizer("OVERWRITE");
1255
1256      Unload() ;
1257      */
1258
1259   }// loader exists
1260   else AliFatal("Loader not available");
1261 }
1262
1263 //__________________________________________________________________
1264 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1265 {
1266   // overloaded method
1267   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1268   if(emcalLoader){
1269
1270     TTree* treeD = emcalLoader->TreeD();
1271     if (!treeD)
1272       {
1273         emcalLoader->MakeDigitsContainer();
1274         treeD = emcalLoader->TreeD();
1275       }
1276
1277     // -- create Digits branch
1278     Int_t bufferSize = 32000;
1279
1280     if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1281       {
1282         triggerBranch->SetAddress(&digits);
1283       }
1284     else
1285       {
1286         treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1287       }
1288
1289     //  treeD->Fill();
1290   }// loader exists
1291   else AliFatal("Loader not available");
1292 }
1293
1294 //__________________________________________________________________
1295 Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1296 {
1297   // Check if cell is defined as dead, so that it is not included
1298   // input is digit
1299   Int_t absId   = digit->GetId();
1300
1301   return IsDead(absId);
1302
1303 }
1304
1305
1306 //__________________________________________________________________
1307 Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1308 {
1309   // Check if cell absID is defined as dead, so that it is not included
1310
1311     AliRunLoader *runLoader = AliRunLoader::Instance();
1312     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1313     if (!emcalLoader)
1314     {
1315       AliFatal("Did not get the  Loader");
1316       return kTRUE;
1317     }
1318
1319     AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1320     if (!caloPed)
1321     {
1322         AliWarning("Could not access pedestal data! No dead channel removal applied");
1323         return kFALSE;
1324     }
1325
1326         // Load Geometry
1327     const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1328     if (!geom)
1329     {
1330       AliFatal("Did not get geometry from EMCALLoader");
1331       return kTRUE; //coverity
1332     }
1333
1334     Int_t iSupMod = -1;
1335     Int_t nModule = -1;
1336     Int_t nIphi   = -1;
1337     Int_t nIeta   = -1;
1338     Int_t iphi    = -1;
1339     Int_t ieta    = -1;
1340
1341     Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1342
1343     if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1344     geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1345
1346     Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1347
1348     if (channelStatus == AliCaloCalibPedestal::kDead)
1349         return kTRUE;
1350     else
1351         return kFALSE;
1352 }
1353