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