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