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