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               if(sdigits->GetEntriesFast() > index[i]){
 479                 AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
 480                 if (tmpdigit){
 481                   curNext = tmpdigit->GetId() ;
 482                 }
 483               }
 484               if(curNext < nextSig) nextSig = curNext ;
 485             }// sdigits exist
 486           } // input loop
 487
 488         }//absID==nextSig
 489
 490         // add the noise now
 491         energy += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
 492         // JLK 26-June-2008
 493         //Now digitize the energy according to the sDigitizer method,
 494         //which merely converts the energy to an integer.  Later we will
 495         //check that the stored value matches our allowed dynamic ranges
 496         digit->SetAmplitude(sDigitizer->Digitize(energy)) ;
 497         AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
 498                          absID, energy, nextSig));
 499       }// Digit pointer exists
 500       else{
 501         Info("Digitizer","Digit pointer is null");
 502       }
 503     } // for(absID = 0; absID < nEMC; absID++)
 504
 505     //ticks->Delete() ;
 506     //delete ticks ;
 507
 508     //JLK is it better to call Clear() here?
 509     delete sdigArray ; //We should not delete its contents
 510
 511     //remove digits below thresholds
 512     // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
 513     // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3
 514     Float_t energy=0;
 515     Float_t timeResolution = 0;
 516     for(i = 0 ; i < nEMC ; i++){
 517       digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
 518       if(digit){
 519         //First get the energy in GeV units.
 520         energy = sDigitizer->Calibrate(digit->GetAmplitude()) ;
 521         //Then digitize using the calibration constants of the ocdb
 522         Float_t ampADC = DigitizeEnergy(energy, digit->GetId())  ;
 523         //if(ampADC>2)printf("Digit energy %f, amp %d, amp cal %d, threshold %d\n",energy,digit->GetAmplitude(),ampADC,fDigitThreshold);
 524         if(ampADC < fDigitThreshold)
 525           digits->RemoveAt(i) ;
 526         else {
 527           timeResolution = GetTimeResolution(energy);
 528           digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
 529         }
 530       }// digit exists
 531     } // digit loop
 532
 533     digits->Compress() ;
 534
 535     Int_t ndigits = digits->GetEntriesFast() ;
 536
 537     //JLK 26-June-2008
 538     //After we have done the summing and digitizing to create the
 539     //digits, now we want to calibrate the resulting amplitude to match
 540     //the dynamic range of our real data.
 541     for (i = 0 ; i < ndigits ; i++) {
 542       digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
 543       if(digit){
 544         digit->SetIndexInList(i) ;
 545         energy = sDigitizer->Calibrate(digit->GetAmplitude()) ;
 546         digit->SetAmplitude(DigitizeEnergy(energy, digit->GetId()) ) ;
 547         //Add delay to time
 548         digit->SetTime(digit->GetTime()+fTimeDelay) ;
 549         // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
 550       }// digit exists
 551     }//Digit loop
 552
 553   }//SDigitizer not null
 554   }
 555   else AliFatal("EMCALLoader is NULL!") ;
 556
 557 }
 558
 559 // //_____________________________________________________________________
 560 Float_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
 561 {
 562   // JLK 26-June-2008
 563   // Returns digitized value of the energy in a cell absId
 564   // using the calibration constants stored in the OCDB
 565   // or default values if no CalibData object is found.
 566   // This effectively converts everything to match the dynamic range
 567   // of the real data we will collect
 568   //
 569   // Load Geometry
 570   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 571
 572   Float_t channel = -999;
 573
 574   if (geom==0){
 575     AliFatal("Did not get geometry from EMCALLoader");
 576   }
 577   else {
 578     Int_t iSupMod = -1;
 579     Int_t nModule  = -1;
 580     Int_t nIphi   = -1;
 581     Int_t nIeta   = -1;
 582     Int_t iphi    = -1;
 583     Int_t ieta    = -1;
 584     Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
 585     if(!bCell)
 586       Error("DigitizeEnergy","Wrong cell id number : AbsId %i ", AbsId) ;
 587     geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 588
 589     if(fCalibData) {
 590       fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 591       fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
 592     }
 593
 594     //channel = static_cast<Int_t> (TMath::Floor( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
 595     channel = (energy + fADCpedestalEC)/fADCchannelEC   ;
 596
 597     if(channel > fNADCEC )
 598       channel =  fNADCEC ;
 599   }// geometry exists
 600
 601   return channel ;
 602
 603 }
 604
 605 //____________________________________________________________________________
 606 void AliEMCALDigitizer::Exec(Option_t *option)
 607 {
 608   // Steering method to process digitization for events
 609   // in the range from fFirstEvent to fLastEvent.
 610   // This range is optionally set by SetEventRange().
 611   // if fLastEvent=-1, then process events until the end.
 612   // by default fLastEvent = fFirstEvent (process only one event)
 613
 614   if (!fInit) { // to prevent overwrite existing file
 615     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
 616     return ;
 617   }
 618
 619   if (strstr(option,"print")) {
 620
 621     Print();
 622     return ;
 623   }
 624
 625   if(strstr(option,"tim"))
 626     gBenchmark->Start("EMCALDigitizer");
 627
 628   AliRunLoader *rl = AliRunLoader::Instance();
 629   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 630   Int_t nEvents = 0;
 631   if(!emcalLoader){
 632     AliFatal("Did not get the  Loader");
 633   }
 634   else{
 635     // Post Digitizer to the white board
 636     emcalLoader->PostDigitizer(this) ;
 637
 638     if (fLastEvent == -1)  {
 639       fLastEvent = rl->GetNumberOfEvents() - 1 ;
 640     }
 641     else if (fManager)
 642       fLastEvent = fFirstEvent ; // what is this ??
 643
 644     nEvents = fLastEvent - fFirstEvent + 1;
 645     Int_t ievent  = -1;
 646
 647     TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit",    32*96);
 648     TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", 32*96);
 649
 650     rl->LoadSDigits("EMCAL");
 651     for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 652
 653       rl->GetEvent(ievent);
 654
 655       Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 656
 657       WriteDigits() ;
 658
 659       //Trigger Digits
 660       //-------------------------------------
 661
 662
 663       Digits2FastOR(digitsTMP, digitsTRG);
 664
 665       WriteDigits(digitsTRG);
 666
 667       (emcalLoader->TreeD())->Fill();
 668
 669       emcalLoader->WriteDigits(   "OVERWRITE");
 670       emcalLoader->WriteDigitizer("OVERWRITE");
 671
 672       Unload();
 673
 674       digitsTRG  ->Delete();
 675       digitsTMP  ->Delete();
 676
 677       //-------------------------------------
 678
 679       if(strstr(option,"deb"))
 680         PrintDigits(option);
 681       if(strstr(option,"table")) gObjectTable->Print();
 682
 683       //increment the total number of Digits per run
 684       fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
 685     }//loop
 686
 687     emcalLoader->CleanDigitizer() ;
 688
 689   }//loader exists
 690
 691   if(strstr(option,"tim")){
 692     gBenchmark->Stop("EMCALDigitizer");
 693     AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event",
 694          gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
 695   }
 696 }
 697
 698 //__________________________________________________________________
 699 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
 700 {
 701   // From F. Blanco
 702   Float_t res = -1;
 703   if (energy > 0) {
 704     res = TMath::Sqrt(fTimeResolutionPar0 +
 705                       fTimeResolutionPar1/(energy*energy) );
 706   }
 707   // parametrization above is for ns. Convert to seconds:
 708   return res*1e-9;
 709 }
 710
 711 //____________________________________________________________________________
 712 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
 713 {
 714   // FEE digits afterburner to produce TRG digits
 715   // we are only interested in the FEE digit deposited energy
 716   // to be converted later into a voltage value
 717
 718   // push the FEE digit to its associated FastOR (numbered from 0:95)
 719   // TRU is in charge of summing module digits
 720
 721   AliRunLoader *runLoader = AliRunLoader::Instance();
 722
 723   AliRun* run = runLoader->GetAliRun();
 724
 725   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
 726   if(!emcalLoader){
 727     AliFatal("Did not get the  Loader");
 728   }
 729   else {
 730     AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(run->GetDetector("EMCAL"));
 731     if(emcal){
 732       AliEMCALGeometry* geom = emcal->GetGeometry();
 733
 734       // build FOR from simulated digits
 735       // and xfer to the corresponding TRU input (mapping)
 736
 737       TClonesArray* digits = emcalLoader->Digits();
 738
 739       TIter NextDigit(digits);
 740       while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
 741       {
 742         Int_t id = digit->GetId();
 743
 744         Int_t trgid;
 745         if (geom && geom->GetFastORIndexFromCellIndex(id , trgid))
 746         {
 747           AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
 748
 749           AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
 750
 751           if (!d)
 752           {
 753             new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
 754             d = (AliEMCALDigit*)digitsTMP->At(trgid);
 755             d->SetId(trgid);
 756           }
 757           else
 758           {
 759             *d = *d + *digit;
 760           }
 761         }
 762       }
 763
 764       if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
 765
 766       Int_t    nSamples = 32;
 767       Int_t *timeSamples = new Int_t[nSamples];
 768
 769       NextDigit = TIter(digitsTMP);
 770       while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
 771       {
 772         if (digit)
 773         {
 774           Int_t     id = digit->GetId();
 775           Float_t time = digit->GetTime();
 776
 777           Double_t depositedEnergy = 0.;
 778           for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
 779
 780           if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
 781
 782           // FIXME: Check digit time!
 783           if (depositedEnergy)
 784           {
 785             DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
 786
 787             for (Int_t j=0;j<nSamples;j++)
 788             {
 789               timeSamples[j] = ((j << 12) & 0xFF000) | (timeSamples[j] & 0xFFF);
 790             }
 791
 792             new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
 793           }
 794         }
 795       }
 796
 797       delete [] timeSamples;
 798     }// AliEMCAL exists
 799     else AliFatal("Could not get AliEMCAL");
 800   }// loader exists
 801
 802 }
 803
 804 //____________________________________________________________________________
 805 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
 806 {
 807         // parameters:
 808         // id: 0..95
 809         const Int_t    reso = 11;      // 11-bit resolution ADC
 810         const Double_t vFSR = 1;       // Full scale input voltage range
 811         const Double_t dNe   = 125;     // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
 812         const Double_t vA   = .136e-6; // CSP output range: 0.136uV/e-
 813         const Double_t rise = 40e-9;   // rise time (10-90%) of the FastOR signal before shaping
 814
 815         const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
 816
 817         Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
 818
 819         TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
 820         signalF.SetParameter( 0,   vV );
 821         signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
 822         signalF.SetParameter( 2, rise );
 823
 824         for (Int_t iTime=0; iTime<nSamples; iTime++)
 825         {
 826                 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
 827                 // probably plan an access to OCDB
 828
 829                 timeSamples[iTime] = int((TMath::Power(2, reso) / vFSR) * signalF.Eval(iTime * kTimeBinWidth) + 0.5);
 830         }
 831 }
 832
 833
 834 //____________________________________________________________________________
 835 //Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks)
 836 //{
 837 //  //  Returns the shortest time among all time ticks
 838 //
 839 //  ticks->Sort() ; //Sort in accordance with times of ticks
 840 //  TIter it(ticks) ;
 841 //  AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
 842 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 843 //
 844 //  AliEMCALTick * t ;
 845 //  while((t=(AliEMCALTick*) it.Next())){
 846 //    if(t->GetTime() < time)  //This tick starts before crossing
 847 //      *ctick+=*t ;
 848 //    else
 849 //      return time ;
 850 //
 851 //    time = ctick->CrossingTime(fTimeThreshold) ;
 852 //  }
 853 //  return time ;
 854 //}
 855 //
 856
 857 //____________________________________________________________________________
 858 Bool_t AliEMCALDigitizer::Init()
 859 {
 860   // Makes all memory allocations
 861   fInit = kTRUE ;
 862   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 863
 864   if ( emcalLoader == 0 ) {
 865     Fatal("Init", "Could not obtain the AliEMCALLoader");
 866     return kFALSE;
 867   }
 868
 869   fFirstEvent = 0 ;
 870   fLastEvent = fFirstEvent ;
 871
 872   if(fManager)
 873     fInput = fManager->GetNinputs() ;
 874   else
 875     fInput           = 1 ;
 876
 877   fInputFileNames  = new TString[fInput] ;
 878   fEventNames      = new TString[fInput] ;
 879   fInputFileNames[0] = GetTitle() ;
 880   fEventNames[0]     = fEventFolderName.Data() ;
 881   Int_t index ;
 882   for (index = 1 ; index < fInput ; index++) {
 883     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 884     TString tempo = fManager->GetInputFolderName(index) ;
 885     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager
 886   }
 887
 888   //to prevent cleaning of this object while GetEvent is called
 889   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 890
 891   //Calibration instance
 892   fCalibData = emcalLoader->CalibData();
 893   return fInit ;
 894 }
 895
 896 //____________________________________________________________________________
 897 void AliEMCALDigitizer::InitParameters()
 898 {
 899   // Parameter initialization for digitizer
 900
 901   // Get the parameters from the OCDB via the loader
 902   AliRunLoader *rl = AliRunLoader::Instance();
 903   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 904   AliEMCALSimParam * simParam = 0x0;
 905   if(emcalLoader) simParam = emcalLoader->SimulationParameters();
 906
 907   if(!simParam){
 908           simParam = AliEMCALSimParam::GetInstance();
 909           AliWarning("Simulation Parameters not available in OCDB?");
 910   }
 911
 912   fMeanPhotonElectron = simParam->GetMeanPhotonElectron();//4400;  // electrons per GeV
 913   fPinNoise           = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
 914   if (fPinNoise < 0.0001 )
 915     Warning("InitParameters", "No noise added\n") ;
 916   fDigitThreshold     = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
 917   fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
 918   fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
 919   fTimeDelay          = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
 920
 921   // These defaults are normally not used.
 922   // Values are read from calibration database instead
 923   fADCchannelEC       = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
 924   fADCpedestalEC      = 0.0 ;  // GeV
 925
 926   fNADCEC          = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
 927
 928   AliDebug(2,Form("Mean Photon Electron %d, Noise %f, Digit Threshold %d,Time Resolution Par0 %g Par1 %g,NADCEC %d",
 929                 fMeanPhotonElectron,fPinNoise,fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
 930
 931   // Not used anymore, remove?
 932   // fTimeSignalLength   = 1.0e-9 ;
 933   // fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
 934
 935 }
 936
 937 //__________________________________________________________________
 938 void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
 939 {
 940   // Allows to produce digits by superimposing background and signal event.
 941   // It is assumed, that headers file with SIGNAL events is opened in
 942   // the constructor.
 943   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed
 944   // Thus we avoid writing (changing) huge and expensive
 945   // backgound files: all output will be writen into SIGNAL, i.e.
 946   // opened in constructor file.
 947   //
 948   // One can open as many files to mix with as one needs.
 949   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
 950   // can be mixed.
 951
 952   if( strcmp(GetName(), "") == 0 )
 953     Init() ;
 954
 955   if(fManager){
 956     Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
 957     return ;
 958   }
 959   // looking for file which contains AliRun
 960   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
 961     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
 962     return ;
 963   }
 964   // looking for the file which contains SDigits
 965   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 966
 967   if(!emcalLoader){
 968     AliFatal("Did not get the  Loader");
 969   }
 970   else{
 971     TString fileName( emcalLoader->GetSDigitsFileName() ) ;
 972     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 973       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
 974     if ( (gSystem->AccessPathName(fileName)) ) {
 975       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
 976       return ;
 977     }
 978     // need to increase the arrays
 979     // MvL: This code only works when fInput == 1, I think.
 980     TString tempo = fInputFileNames[fInput-1] ;
 981     delete [] fInputFileNames ;
 982     fInputFileNames = new TString[fInput+1] ;
 983     fInputFileNames[fInput-1] = tempo ;
 984
 985     tempo = fEventNames[fInput-1] ;
 986     delete [] fEventNames ;
 987     fEventNames = new TString[fInput+1] ;
 988     fEventNames[fInput-1] = tempo ;
 989
 990     fInputFileNames[fInput] = alirunFileName ;
 991     fEventNames[fInput]     = eventFolderName ;
 992     fInput++ ;
 993
 994   }// Loader exits
 995
 996 }
 997
 998 //__________________________________________________________________
 999 void AliEMCALDigitizer::Print1(Option_t * option)
1000 { // 19-nov-04 - just for convinience
1001   Print();
1002   PrintDigits(option);
1003 }
1004
1005 //__________________________________________________________________
1006 void AliEMCALDigitizer::Print(Option_t*)const
1007 {
1008   // Print Digitizer's parameters
1009   printf("Print: \n------------------- %s -------------", GetName() ) ;
1010   if( strcmp(fEventFolderName.Data(), "") != 0 ){
1011     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
1012
1013     Int_t nStreams ;
1014     if (fManager)
1015       nStreams =  GetNInputStreams() ;
1016     else
1017       nStreams = fInput ;
1018
1019     AliRunLoader *rl=0;
1020
1021     Int_t index = 0 ;
1022     for (index = 0 ; index < nStreams ; index++) {
1023       TString tempo(fEventNames[index]) ;
1024       tempo += index ;
1025       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1026         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1027       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1028       if(emcalLoader){
1029         TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1030         if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1031           fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
1032         printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1033       }//loader
1034     }
1035
1036     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1037
1038     if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1039     else printf("\nNULL LOADER");
1040
1041     printf("\nWith following parameters:\n") ;
1042
1043     printf("    Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise) ;
1044     printf("    Threshold  in Tower  (fDigitThreshold) = %d\n", fDigitThreshold)  ;
1045     printf("---------------------------------------------------\n")  ;
1046   }
1047   else
1048     printf("Print: AliEMCALDigitizer not initialized") ;
1049 }
1050
1051 //__________________________________________________________________
1052 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1053 {
1054   //utility method for printing digit information
1055
1056   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1057   if(emcalLoader){
1058     TClonesArray * digits  = emcalLoader->Digits() ;
1059     TClonesArray * sdigits = emcalLoader->SDigits() ;
1060
1061     printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1062     printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1063
1064     if(strstr(option,"all")){
1065       //loop over digits
1066       AliEMCALDigit * digit;
1067       printf("\nEMC digits (with primaries):\n")  ;
1068       printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
1069       Int_t index ;
1070       for (index = 0 ; index < digits->GetEntries() ; index++) {
1071         digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1072         if(digit){
1073           printf("\n%6d  %8f    %6.5e %4d      %2d : ",
1074                  digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1075           Int_t iprimary;
1076           for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1077             printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1078           }
1079         }// digit exists
1080       }// loop
1081     }
1082     printf("\n");
1083   }// loader exists
1084   else printf("NULL LOADER, cannot print\n");
1085 }
1086
1087 //__________________________________________________________________
1088 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1089 {
1090   // Calculates the time signal generated by noise
1091   //PH  Info("TimeOfNoise", "Change me") ;
1092   return gRandom->Rndm() * 1.28E-5;
1093 }
1094
1095 //__________________________________________________________________
1096 void AliEMCALDigitizer::Unload()
1097 {
1098   // Unloads the SDigits and Digits
1099   AliRunLoader *rl=0;
1100
1101   Int_t i ;
1102   for(i = 1 ; i < fInput ; i++){
1103     TString tempo(fEventNames[i]) ;
1104     tempo += i;
1105     if ((rl = AliRunLoader::GetRunLoader(tempo)))
1106       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1107   }
1108   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1109   if(emcalLoader)emcalLoader->UnloadDigits() ;
1110   else AliFatal("Did not get the loader");
1111 }
1112
1113 //_________________________________________________________________________________________
1114 void AliEMCALDigitizer::WriteDigits()
1115 { // Makes TreeD in the output file.
1116   // Check if branch already exists:
1117   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1118   //      already existing branches.
1119   //   else creates branch with Digits, named "EMCAL", title "...",
1120   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1121   //      and names of files, from which digits are made.
1122
1123   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1124
1125   if(emcalLoader){
1126     const TClonesArray * digits = emcalLoader->Digits() ;
1127     TTree * treeD = emcalLoader->TreeD();
1128     if ( !treeD ) {
1129       emcalLoader->MakeDigitsContainer();
1130       treeD = emcalLoader->TreeD();
1131     }
1132
1133     // -- create Digits branch
1134     Int_t bufferSize = 32000 ;
1135     TBranch * digitsBranch = 0;
1136     if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1137       digitsBranch->SetAddress(&digits);
1138       AliWarning("Digits Branch already exists. Not all digits will be visible");
1139     }
1140     else
1141       treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1142     //digitsBranch->SetTitle(fEventFolderName);
1143
1144     //  treeD->Fill() ;
1145     /*
1146      emcalLoader->WriteDigits("OVERWRITE");
1147      emcalLoader->WriteDigitizer("OVERWRITE");
1148
1149      Unload() ;
1150      */
1151
1152   }// loader exists
1153   else AliFatal("Loader not available");
1154 }
1155
1156 //__________________________________________________________________
1157 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1158 { // overloaded method
1159         AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1160         if(emcalLoader){
1161
1162     TTree* treeD = emcalLoader->TreeD();
1163     if (!treeD)
1164     {
1165       emcalLoader->MakeDigitsContainer();
1166       treeD = emcalLoader->TreeD();
1167     }
1168
1169     // -- create Digits branch
1170     Int_t bufferSize = 32000;
1171
1172     if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1173     {
1174       triggerBranch->SetAddress(&digits);
1175     }
1176     else
1177     {
1178       treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1179     }
1180
1181     //  treeD->Fill();
1182   }// loader exists
1183   else AliFatal("Loader not available");
1184 }
1185
1186 //__________________________________________________________________
1187 void AliEMCALDigitizer::Browse(TBrowser* b)
1188 {
1189   TTask::Browse(b);
1190 }