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