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