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