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 <cassert>
  70
  71 // --- AliRoot header files ---
  72 #include "AliLog.h"
  73 #include "AliRun.h"
  74 #include "AliRunDigitizer.h"
  75 #include "AliRunLoader.h"
  76 #include "AliCDBManager.h"
  77 #include "AliCDBEntry.h"
  78 #include "AliEMCALDigit.h"
  79 #include "AliEMCAL.h"
  80 #include "AliEMCALLoader.h"
  81 #include "AliEMCALDigitizer.h"
  82 #include "AliEMCALSDigitizer.h"
  83 #include "AliEMCALGeometry.h"
  84 #include "AliEMCALTick.h"
  85
  86 ClassImp(AliEMCALDigitizer)
  87
  88
  89 //____________________________________________________________________________
  90 AliEMCALDigitizer::AliEMCALDigitizer()
  91   : AliDigitizer("",""),
  92     fDefaultInit(kTRUE),
  93     fDigitsInRun(0),
  94     fInit(0),
  95     fInput(0),
  96     fInputFileNames(0x0),
  97     fEventNames(0x0),
  98     fDigitThreshold(0),
  99     fMeanPhotonElectron(0),
 100     fPedestal(0),
 101     fSlope(0),
 102     fPinNoise(0),
 103     fTimeResolution(0),
 104     fTimeThreshold(0),
 105     fTimeSignalLength(0),
 106     fADCchannelEC(0),
 107     fADCpedestalEC(0),
 108     fNADCEC(0),
 109     fEventFolderName(""),
 110     fFirstEvent(0),
 111     fLastEvent(0),
 112     fCalibData(0x0)
 113 {
 114   // ctor
 115   InitParameters() ;
 116   fManager = 0 ;                     // We work in the standalone mode
 117 }
 118
 119 //____________________________________________________________________________
 120 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
 121   : AliDigitizer("EMCAL"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
 122     fDefaultInit(kFALSE),
 123     fDigitsInRun(0),
 124     fInit(0),
 125     fInput(0),
 126     fInputFileNames(0),
 127     fEventNames(0),
 128     fDigitThreshold(0),
 129     fMeanPhotonElectron(0),
 130     fPedestal(0),
 131     fSlope(0),
 132     fPinNoise(0),
 133     fTimeResolution(0),
 134     fTimeThreshold(0),
 135     fTimeSignalLength(0),
 136     fADCchannelEC(0),
 137     fADCpedestalEC(0),
 138     fNADCEC(0),
 139     fEventFolderName(eventFolderName),
 140     fFirstEvent(0),
 141     fLastEvent(0),
 142     fCalibData(0x0)
 143 {
 144   // ctor
 145   InitParameters() ;
 146   Init() ;
 147   fManager = 0 ;                     // We work in the standalone mode
 148 }
 149
 150 //____________________________________________________________________________
 151 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
 152   : AliDigitizer(d.GetName(),d.GetTitle()),
 153     fDefaultInit(d.fDefaultInit),
 154     fDigitsInRun(d.fDigitsInRun),
 155     fInit(d.fInit),
 156     fInput(d.fInput),
 157     fInputFileNames(d.fInputFileNames),
 158     fEventNames(d.fEventNames),
 159     fDigitThreshold(d.fDigitThreshold),
 160     fMeanPhotonElectron(d.fMeanPhotonElectron),
 161     fPedestal(d.fPedestal),
 162     fSlope(d.fSlope),
 163     fPinNoise(d.fPinNoise),
 164     fTimeResolution(d.fTimeResolution),
 165     fTimeThreshold(d.fTimeThreshold),
 166     fTimeSignalLength(d.fTimeSignalLength),
 167     fADCchannelEC(d.fADCchannelEC),
 168     fADCpedestalEC(d.fADCpedestalEC),
 169     fNADCEC(d.fNADCEC),
 170     fEventFolderName(d.fEventFolderName),
 171     fFirstEvent(d.fFirstEvent),
 172     fLastEvent(d.fLastEvent),
 173     fCalibData(d.fCalibData)
 174 {
 175   // copyy ctor
 176  }
 177
 178 //____________________________________________________________________________
 179 AliEMCALDigitizer::AliEMCALDigitizer(AliRunDigitizer * rd)
 180   : AliDigitizer(rd,"EMCAL"+AliConfig::Instance()->GetDigitizerTaskName()),
 181     fDefaultInit(kFALSE),
 182     fDigitsInRun(0),
 183     fInit(0),
 184     fInput(0),
 185     fInputFileNames(0),
 186     fEventNames(0),
 187     fDigitThreshold(0.),
 188     fMeanPhotonElectron(0),
 189     fPedestal(0),
 190     fSlope(0.),
 191     fPinNoise(0),
 192     fTimeResolution(0.),
 193     fTimeThreshold(0),
 194     fTimeSignalLength(0),
 195     fADCchannelEC(0),
 196     fADCpedestalEC(0),
 197     fNADCEC(0),
 198     fEventFolderName(0),
 199     fFirstEvent(0),
 200     fLastEvent(0),
 201     fCalibData(0x0)
 202 {
 203   // ctor Init() is called by RunDigitizer
 204   fManager = rd ;
 205   fEventFolderName = fManager->GetInputFolderName(0) ;
 206   SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
 207   InitParameters() ;
 208 }
 209
 210 //____________________________________________________________________________
 211   AliEMCALDigitizer::~AliEMCALDigitizer()
 212 {
 213   //dtor
 214   if (AliRunLoader::Instance()) {
 215     AliLoader *emcalLoader=0;
 216     if ((emcalLoader = AliRunLoader::Instance()->GetDetectorLoader("EMCAL")))
 217       emcalLoader->CleanDigitizer();
 218   }
 219   else
 220     AliDebug(1," no runloader present");
 221   delete [] fInputFileNames ;
 222   delete [] fEventNames ;
 223
 224 }
 225
 226 //____________________________________________________________________________
 227 void AliEMCALDigitizer::Digitize(Int_t event)
 228 {
 229
 230   // Makes the digitization of the collected summable digits
 231   // for this it first creates the array of all EMCAL modules
 232   // filled with noise and after that adds contributions from
 233   // SDigits. This design helps to avoid scanning over the
 234   // list of digits to add  contribution of any new SDigit.
 235   //
 236   // JLK 26-Jun-2008
 237   // Note that SDigit energy info is stored as an amplitude, so we
 238   // must call the Calibrate() method of the SDigitizer to convert it
 239   // back to an energy in GeV before adding it to the Digit
 240   //
 241   static int nEMC=0; //max number of digits possible
 242
 243   AliRunLoader *rl = AliRunLoader::Instance();
 244   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 245   Int_t readEvent = event ;
 246   // fManager is data member from AliDigitizer
 247   if (fManager)
 248     readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
 249   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 250                   readEvent, GetTitle(), fEventFolderName.Data())) ;
 251   rl->GetEvent(readEvent);
 252
 253   TClonesArray * digits = emcalLoader->Digits() ;
 254   digits->Delete() ;  //correct way to clear array when memory is
 255                       //allocated by objects stored in array
 256
 257   // Load Geometry
 258   AliEMCALGeometry *geom = 0;
 259   if (rl->GetAliRun()) {
 260     AliEMCAL * emcal  = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
 261     geom = emcal->GetGeometry();
 262   }
 263   else
 264     AliFatal("Could not get AliRun from runLoader");
 265
 266   nEMC = geom->GetNCells();
 267   AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
 268
 269   Int_t absID ;
 270
 271   digits->Expand(nEMC) ;
 272
 273   // get first the sdigitizer from the tasks list (must have same name as the digitizer)
 274   if ( !emcalLoader->SDigitizer() )
 275     emcalLoader->LoadSDigitizer();
 276   AliEMCALSDigitizer * sDigitizer = dynamic_cast<AliEMCALSDigitizer*>(emcalLoader->SDigitizer());
 277
 278   if ( !sDigitizer )
 279     Fatal("Digitize", "SDigitizer with name %s %s not found", fEventFolderName.Data(), GetTitle() ) ;
 280
 281   //take all the inputs to add together and load the SDigits
 282   TObjArray * sdigArray = new TObjArray(fInput) ;
 283   sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
 284   Int_t i ;
 285
 286   for(i = 1 ; i < fInput ; i++){
 287     TString tempo(fEventNames[i]) ;
 288     tempo += i ;
 289
 290     AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ;
 291
 292     if (rl2==0)
 293       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
 294
 295     if (fManager)
 296       readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
 297     Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
 298     rl2->LoadSDigits();
 299     rl2->GetEvent(readEvent);
 300     AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
 301     sdigArray->AddAt(emcalLoader2->SDigits(), i) ;
 302   }
 303
 304   //Find the first tower with signal
 305   Int_t nextSig = nEMC + 1 ;
 306   TClonesArray * sdigits ;
 307   for(i = 0 ; i < fInput ; i++){
 308     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 309     if ( !sdigits->GetEntriesFast() )
 310       continue ;
 311     Int_t curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(0))->GetId() ;
 312      if(curNext < nextSig)
 313        nextSig = curNext ;
 314      AliDebug(1,Form("input %i : #sdigits %i \n",
 315                      i, sdigits->GetEntriesFast()));
 316   }
 317   AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
 318
 319   TArrayI index(fInput) ;
 320   index.Reset() ;  //Set all indexes to zero
 321
 322   AliEMCALDigit * digit ;
 323   AliEMCALDigit * curSDigit ;
 324
 325   //  TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ;
 326
 327   //Put Noise contribution
 328   for(absID = 0; absID < nEMC; absID++){ // Nov 30, 2006 by PAI; was from 1 to nEMC
 329     Float_t energy = 0 ;
 330     // amplitude set to zero, noise will be added later
 331     new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0, TimeOfNoise() ); // absID-1->absID
 332     //look if we have to add signal?
 333     digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
 334
 335     if(absID==nextSig){
 336       //Add SDigits from all inputs
 337       //      ticks->Clear() ;
 338       //Int_t contrib = 0 ;
 339
 340       //Follow PHOS and comment out this timing model til a better one
 341       //can be developed - JLK 28-Apr-2008
 342
 343       //Float_t a = digit->GetAmp() ;
 344       //Float_t b = TMath::Abs( a /fTimeSignalLength) ;
 345       //Mark the beginning of the signal
 346       //new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime(),0, b);
 347       //Mark the end of the signal
 348       //new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime()+fTimeSignalLength, -a, -b);
 349
 350       // Calculate time as time of the largest digit
 351       Float_t time = digit->GetTime() ;
 352       Float_t aTime= digit->GetAmp() ;
 353
 354       // loop over input
 355       for(i = 0; i< fInput ; i++){  //loop over (possible) merge sources
 356         if(dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 357           curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 358         else
 359           curSDigit = 0 ;
 360         //May be several digits will contribute from the same input
 361         while(curSDigit && (curSDigit->GetId() == absID)){
 362           //Shift primary to separate primaries belonging different inputs
 363           Int_t primaryoffset ;
 364           if(fManager)
 365             primaryoffset = fManager->GetMask(i) ;
 366           else
 367             primaryoffset = i ;
 368           curSDigit->ShiftPrimary(primaryoffset) ;
 369
 370           //Remove old timing model - JLK 28-April-2008
 371           //a = curSDigit->GetAmp() ;
 372           //b = a /fTimeSignalLength ;
 373           //new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime(),0, b);
 374           //new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
 375           if(curSDigit->GetAmp()>aTime) {
 376             aTime = curSDigit->GetAmp();
 377             time = curSDigit->GetTime();
 378           }
 379
 380           *digit = *digit + *curSDigit ;  //adds amplitudes of each digit
 381
 382           index[i]++ ;
 383           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 384             curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 385           else
 386             curSDigit = 0 ;
 387         }
 388       }
 389       //Here we convert the summed amplitude to an energy in GeV
 390       energy = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
 391       // add fluctuations for photo-electron creation
 392       energy *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
 393
 394       //calculate and set time
 395       //New timing model needed - JLK 28-April-2008
 396       //Float_t time = FrontEdgeTime(ticks) ;
 397       digit->SetTime(time) ;
 398
 399       //Find next signal module
 400       nextSig = nEMC + 1 ;
 401       for(i = 0 ; i < fInput ; i++){
 402         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 403         Int_t curNext = nextSig ;
 404         if(sdigits->GetEntriesFast() > index[i] ){
 405           curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]))->GetId() ;
 406         }
 407         if(curNext < nextSig) nextSig = curNext ;
 408       }
 409     }
 410     // add the noise now
 411     energy += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
 412     // JLK 26-June-2008
 413     //Now digitize the energy according to the sDigitizer method,
 414     //which merely converts the energy to an integer.  Later we will
 415     //check that the stored value matches our allowed dynamic ranges
 416     digit->SetAmp(sDigitizer->Digitize(energy)) ;
 417     AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
 418                      absID, energy, nextSig));
 419   } // for(absID = 0; absID < nEMC; absID++)
 420
 421   //ticks->Delete() ;
 422   //delete ticks ;
 423
 424   //JLK is it better to call Clear() here?
 425   delete sdigArray ; //We should not delete its contents
 426
 427   //remove digits below thresholds
 428   for(i = 0 ; i < nEMC ; i++){
 429     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
 430     Float_t threshold = fDigitThreshold ; //this is in GeV
 431     //need to calibrate digit amplitude to energy in GeV for comparison
 432     if(sDigitizer->Calibrate( digit->GetAmp() ) < threshold)
 433       digits->RemoveAt(i) ;
 434     else
 435       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
 436   }
 437
 438   digits->Compress() ;
 439
 440   Int_t ndigits = digits->GetEntriesFast() ;
 441
 442   //JLK 26-June-2008
 443   //After we have done the summing and digitizing to create the
 444   //digits, now we want to calibrate the resulting amplitude to match
 445   //the dynamic range of our real data.
 446   Float_t energy=0;
 447   for (i = 0 ; i < ndigits ; i++) {
 448     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
 449     digit->SetIndexInList(i) ;
 450     energy = sDigitizer->Calibrate(digit->GetAmp()) ;
 451     digit->SetAmp(DigitizeEnergy(energy, digit->GetId()) ) ;
 452   }
 453
 454 }
 455
 456 // //_____________________________________________________________________
 457 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
 458 {
 459   // JLK 26-June-2008
 460   // Returns digitized value of the energy in a cell absId
 461   // using the calibration constants stored in the OCDB
 462   // or default values if no CalibData object is found.
 463   // This effectively converts everything to match the dynamic range
 464   // of the real data we will collect
 465   //
 466   // Load Geometry
 467   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 468
 469   if (geom==0)
 470     AliFatal("Did not get geometry from EMCALLoader");
 471
 472   Int_t iSupMod = -1;
 473   Int_t nModule  = -1;
 474   Int_t nIphi   = -1;
 475   Int_t nIeta   = -1;
 476   Int_t iphi    = -1;
 477   Int_t ieta    = -1;
 478   Int_t channel = -999;
 479
 480   Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
 481   if(!bCell)
 482     Error("DigitizeEnergy","Wrong cell id number : AbsId %i ", AbsId) ;
 483   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 484
 485   if(fCalibData) {
 486     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 487     fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
 488   }
 489
 490   channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
 491
 492   if(channel > fNADCEC )
 493     channel =  fNADCEC ;
 494   return channel ;
 495
 496 }
 497
 498 //____________________________________________________________________________
 499 void AliEMCALDigitizer::Exec(Option_t *option)
 500 {
 501   // Steering method to process digitization for events
 502   // in the range from fFirstEvent to fLastEvent.
 503   // This range is optionally set by SetEventRange().
 504   // if fLastEvent=-1, then process events until the end.
 505   // by default fLastEvent = fFirstEvent (process only one event)
 506
 507   if (!fInit) { // to prevent overwrite existing file
 508     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
 509     return ;
 510   }
 511
 512   if (strstr(option,"print")) {
 513
 514     Print();
 515     return ;
 516   }
 517
 518   if(strstr(option,"tim"))
 519     gBenchmark->Start("EMCALDigitizer");
 520
 521   AliRunLoader *rl = AliRunLoader::Instance();
 522   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 523
 524   // Post Digitizer to the white board
 525   emcalLoader->PostDigitizer(this) ;
 526
 527   if (fLastEvent == -1)  {
 528     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 529   }
 530   else if (fManager)
 531     fLastEvent = fFirstEvent ; // what is this ??
 532
 533   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 534   Int_t ievent;
 535
 536   rl->LoadSDigits("EMCAL");
 537   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 538
 539     rl->GetEvent(ievent);
 540
 541     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 542
 543     WriteDigits() ;
 544
 545     if(strstr(option,"deb"))
 546       PrintDigits(option);
 547     if(strstr(option,"table")) gObjectTable->Print();
 548
 549     //increment the total number of Digits per run
 550     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
 551   }
 552
 553   emcalLoader->CleanDigitizer() ;
 554
 555   if(strstr(option,"tim")){
 556     gBenchmark->Stop("EMCALDigitizer");
 557     AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event",
 558          gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
 559   }
 560 }
 561
 562 //____________________________________________________________________________
 563 Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks)
 564 {
 565   //  Returns the shortest time among all time ticks
 566
 567   ticks->Sort() ; //Sort in accordance with times of ticks
 568   TIter it(ticks) ;
 569   AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
 570   Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 571
 572   AliEMCALTick * t ;
 573   while((t=(AliEMCALTick*) it.Next())){
 574     if(t->GetTime() < time)  //This tick starts before crossing
 575       *ctick+=*t ;
 576     else
 577       return time ;
 578
 579     time = ctick->CrossingTime(fTimeThreshold) ;
 580   }
 581   return time ;
 582 }
 583
 584 //____________________________________________________________________________
 585 Bool_t AliEMCALDigitizer::Init()
 586 {
 587   // Makes all memory allocations
 588   fInit = kTRUE ;
 589   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 590
 591   if ( emcalLoader == 0 ) {
 592     Fatal("Init", "Could not obtain the AliEMCALLoader");
 593     return kFALSE;
 594   }
 595
 596   fFirstEvent = 0 ;
 597   fLastEvent = fFirstEvent ;
 598
 599   if(fManager)
 600     fInput = fManager->GetNinputs() ;
 601   else
 602     fInput           = 1 ;
 603
 604   fInputFileNames  = new TString[fInput] ;
 605   fEventNames      = new TString[fInput] ;
 606   fInputFileNames[0] = GetTitle() ;
 607   fEventNames[0]     = fEventFolderName.Data() ;
 608   Int_t index ;
 609   for (index = 1 ; index < fInput ; index++) {
 610     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 611     TString tempo = fManager->GetInputFolderName(index) ;
 612     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager
 613   }
 614
 615   //to prevent cleaning of this object while GetEvent is called
 616   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 617
 618   //Calibration instance
 619   fCalibData = emcalLoader->CalibData();
 620   return fInit ;
 621 }
 622
 623 //____________________________________________________________________________
 624 void AliEMCALDigitizer::InitParameters()
 625 {
 626   // Parameter initialization for digitizer
 627   // Tune parameters - 24-nov-04; Apr 29, 2007
 628   // New parameters JLK 14-Apr-2008
 629
 630   fMeanPhotonElectron = 4400;  // electrons per GeV
 631   fPinNoise           = 0.037; // pin noise in GEV from analysis test beam data
 632   if (fPinNoise == 0. )
 633     Warning("InitParameters", "No noise added\n") ;
 634   fDigitThreshold     = fPinNoise * 3; // 3 * sigma
 635   fTimeResolution     = 0.3e-9 ; // 300 psc
 636   fTimeSignalLength   = 1.0e-9 ;
 637
 638   // These defaults are normally not used.
 639   // Values are read from calibration database instead
 640   fADCchannelEC    = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
 641   fADCpedestalEC   = 0.0 ;  // GeV
 642   fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
 643
 644   fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
 645
 646 }
 647
 648 //__________________________________________________________________
 649 void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
 650 {
 651   // Allows to produce digits by superimposing background and signal event.
 652   // It is assumed, that headers file with SIGNAL events is opened in
 653   // the constructor.
 654   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed
 655   // Thus we avoid writing (changing) huge and expensive
 656   // backgound files: all output will be writen into SIGNAL, i.e.
 657   // opened in constructor file.
 658   //
 659   // One can open as many files to mix with as one needs.
 660   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
 661   // can be mixed.
 662
 663   if( strcmp(GetName(), "") == 0 )
 664     Init() ;
 665
 666   if(fManager){
 667     Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
 668     return ;
 669   }
 670   // looking for file which contains AliRun
 671   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
 672     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
 673     return ;
 674   }
 675   // looking for the file which contains SDigits
 676   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 677   TString fileName( emcalLoader->GetSDigitsFileName() ) ;
 678     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 679       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
 680     if ( (gSystem->AccessPathName(fileName)) ) {
 681       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
 682       return ;
 683     }
 684     // need to increase the arrays
 685     // MvL: This code only works when fInput == 1, I think.
 686     TString tempo = fInputFileNames[fInput-1] ;
 687     delete [] fInputFileNames ;
 688     fInputFileNames = new TString[fInput+1] ;
 689     fInputFileNames[fInput-1] = tempo ;
 690
 691     tempo = fEventNames[fInput-1] ;
 692     delete [] fEventNames ;
 693     fEventNames = new TString[fInput+1] ;
 694     fEventNames[fInput-1] = tempo ;
 695
 696     fInputFileNames[fInput] = alirunFileName ;
 697     fEventNames[fInput]     = eventFolderName ;
 698     fInput++ ;
 699 }
 700
 701 //__________________________________________________________________
 702 void AliEMCALDigitizer::Print1(Option_t * option)
 703 { // 19-nov-04 - just for convinience
 704   Print();
 705   PrintDigits(option);
 706 }
 707
 708 //__________________________________________________________________
 709 void AliEMCALDigitizer::Print(Option_t*)const
 710 {
 711   // Print Digitizer's parameters
 712   printf("Print: \n------------------- %s -------------", GetName() ) ;
 713   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 714     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 715
 716     Int_t nStreams ;
 717     if (fManager)
 718       nStreams =  GetNInputStreams() ;
 719     else
 720       nStreams = fInput ;
 721
 722     Int_t index = 0 ;
 723
 724     AliRunLoader *rl=0;
 725
 726     for (index = 0 ; index < nStreams ; index++) {
 727       TString tempo(fEventNames[index]) ;
 728       tempo += index ;
 729       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
 730         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
 731       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 732       TString fileName( emcalLoader->GetSDigitsFileName() ) ;
 733       if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 734         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
 735       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
 736     }
 737
 738     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 739
 740     printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
 741
 742     printf("\nWith following parameters:\n") ;
 743
 744     printf("    Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise) ;
 745     printf("    Threshold  in EMC  (fDigitThreshold) = %f\n", fDigitThreshold)  ;
 746     printf("---------------------------------------------------\n")  ;
 747   }
 748   else
 749     printf("Print: AliEMCALDigitizer not initialized") ;
 750 }
 751
 752 //__________________________________________________________________
 753 void AliEMCALDigitizer::PrintDigits(Option_t * option)
 754 {
 755   //utility method for printing digit information
 756
 757   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 758   TClonesArray * digits  = emcalLoader->Digits() ;
 759   TClonesArray * sdigits = emcalLoader->SDigits() ;
 760
 761   printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
 762   printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
 763
 764   if(strstr(option,"all")){
 765     //loop over digits
 766     AliEMCALDigit * digit;
 767     printf("\nEMC digits (with primaries):\n")  ;
 768     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 769     Int_t index ;
 770     for (index = 0 ; index < digits->GetEntries() ; index++) {
 771       digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
 772       printf("\n%6d  %8d    %6.5e %4d      %2d : ",
 773               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 774       Int_t iprimary;
 775       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 776         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 777       }
 778     }
 779   }
 780   printf("\n");
 781 }
 782
 783 //__________________________________________________________________
 784 Float_t AliEMCALDigitizer::TimeOfNoise(void)
 785 {
 786   // Calculates the time signal generated by noise
 787   //PH  Info("TimeOfNoise", "Change me") ;
 788   return gRandom->Rndm() * 1.28E-5;
 789 }
 790
 791 //__________________________________________________________________
 792 void AliEMCALDigitizer::Unload()
 793 {
 794   // Unloads the SDigits and Digits
 795   AliRunLoader *rl=0;
 796
 797   Int_t i ;
 798   for(i = 1 ; i < fInput ; i++){
 799     TString tempo(fEventNames[i]) ;
 800     tempo += i;
 801     if ((rl = AliRunLoader::GetRunLoader(tempo)))
 802       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
 803   }
 804   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 805   emcalLoader->UnloadDigits() ;
 806 }
 807
 808 //_________________________________________________________________________________________
 809 void AliEMCALDigitizer::WriteDigits()
 810 {
 811
 812   // Makes TreeD in the output file.
 813   // Check if branch already exists:
 814   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
 815   //      already existing branches.
 816   //   else creates branch with Digits, named "EMCAL", title "...",
 817   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
 818   //      and names of files, from which digits are made.
 819
 820   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 821
 822   const TClonesArray * digits = emcalLoader->Digits() ;
 823   TTree * treeD = emcalLoader->TreeD();
 824   if ( !treeD ) {
 825     emcalLoader->MakeDigitsContainer();
 826     treeD = emcalLoader->TreeD();
 827   }
 828
 829   // -- create Digits branch
 830   Int_t bufferSize = 32000 ;
 831   TBranch * digitsBranch = 0;
 832   if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
 833     digitsBranch->SetAddress(&digits);
 834     AliWarning("Digits Branch already exists. Not all digits will be visible");
 835   }
 836   else
 837     treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
 838   //digitsBranch->SetTitle(fEventFolderName);
 839   treeD->Fill() ;
 840
 841   emcalLoader->WriteDigits("OVERWRITE");
 842   emcalLoader->WriteDigitizer("OVERWRITE");
 843
 844   Unload() ;
 845
 846 }
 847
 848 //__________________________________________________________________
 849 void AliEMCALDigitizer::Browse(TBrowser* b)
 850 {
 851   TTask::Browse(b);
 852 }