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