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() ;  //JLK why is this created then deleted?
 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
 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 eTime= 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()>eTime) {
 378             eTime = curSDigit->GetAmp();
 379             time = curSDigit->GetTime();
 380           }
 381
 382           *digit = *digit + *curSDigit ;  //add energies
 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       // add fluctuations for photo-electron creation
 392       amp = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
 393       amp *= 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     amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
 413     digit->SetAmp(sDigitizer->Digitize(amp)) ;
 414     AliDebug(10,Form(" absID %5i amp %f nextSig %5i\n",
 415                      absID, amp, nextSig));
 416   } // for(absID = 1; absID <= nEMC; absID++)
 417
 418   //ticks->Delete() ;
 419   //delete ticks ;
 420
 421   delete sdigArray ; //We should not delete its contents
 422
 423   //remove digits below thresholds
 424   for(i = 0 ; i < nEMC ; i++){
 425     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
 426     Float_t threshold = fDigitThreshold ;
 427     if(sDigitizer->Calibrate( digit->GetAmp() ) < threshold)
 428       digits->RemoveAt(i) ;
 429     else
 430       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
 431   }
 432
 433   digits->Compress() ;
 434
 435   Int_t ndigits = digits->GetEntriesFast() ;
 436
 437   //Set indexes in list of digits and fill hists.
 438   AliEMCALHistoUtilities::FillH1(fHists, 0, Double_t(ndigits));
 439   Float_t energy=0., esum=0.;
 440   for (i = 0 ; i < ndigits ; i++) {
 441     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
 442     digit->SetIndexInList(i) ;
 443     energy = sDigitizer->Calibrate(digit->GetAmp()) ;
 444     esum += energy;
 445     digit->SetAmp(DigitizeEnergy(energy, digit->GetId()) ) ; // for what ??
 446     AliEMCALHistoUtilities::FillH1(fHists, 2, double(digit->GetAmp()));
 447     AliEMCALHistoUtilities::FillH1(fHists, 3, double(energy));
 448     AliEMCALHistoUtilities::FillH1(fHists, 4, double(digit->GetId()));
 449     //    if(digit->GetId() == nEMC) {
 450     //  printf(" i %i \n", i );
 451     //  digit->Dump();
 452     //  assert(0);
 453     //}
 454   }
 455   AliEMCALHistoUtilities::FillH1(fHists, 1, esum);
 456 }
 457
 458 // //_____________________________________________________________________
 459 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
 460 {
 461   // Returns digitized value of the energy in a cell absId
 462
 463   // Load Geometry
 464   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 465
 466   if (geom==0)
 467     AliFatal("Did not get geometry from EMCALLoader");
 468
 469   Int_t iSupMod = -1;
 470   Int_t nModule  = -1;
 471   Int_t nIphi   = -1;
 472   Int_t nIeta   = -1;
 473   Int_t iphi    = -1;
 474   Int_t ieta    = -1;
 475   Int_t channel = -999;
 476
 477   Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
 478   if(!bCell)
 479     Error("DigitizeEnergy","Wrong cell id number : AbsId %i ", AbsId) ;
 480   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 481
 482   if(fCalibData) {
 483     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 484     fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
 485   }
 486
 487   channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
 488
 489   if(channel > fNADCEC )
 490     channel =  fNADCEC ;
 491   return channel ;
 492
 493 }
 494
 495 //____________________________________________________________________________
 496 void AliEMCALDigitizer::Exec(Option_t *option)
 497 {
 498   // Steering method to process digitization for events
 499   // in the range from fFirstEvent to fLastEvent.
 500   // This range is optionally set by SetEventRange().
 501   // if fLastEvent=-1, then process events until the end.
 502   // by default fLastEvent = fFirstEvent (process only one event)
 503
 504   if (!fInit) { // to prevent overwrite existing file
 505     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
 506     return ;
 507   }
 508
 509   if (strstr(option,"print")) {
 510
 511     Print();
 512     return ;
 513   }
 514
 515   if(strstr(option,"tim"))
 516     gBenchmark->Start("EMCALDigitizer");
 517
 518   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 519   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 520
 521   // Post Digitizer to the white board
 522   emcalLoader->PostDigitizer(this) ;
 523
 524   if (fLastEvent == -1)  {
 525     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 526   }
 527   else if (fManager)
 528     fLastEvent = fFirstEvent ; // what is this ??
 529
 530   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 531   Int_t ievent;
 532
 533   rl->LoadSDigits("EMCAL");
 534   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 535
 536     rl->GetEvent(ievent);
 537
 538     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 539
 540     WriteDigits() ;
 541
 542     if(strstr(option,"deb"))
 543       PrintDigits(option);
 544     if(strstr(option,"table")) gObjectTable->Print();
 545
 546     //increment the total number of Digits per run
 547     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
 548   }
 549
 550   emcalLoader->CleanDigitizer() ;
 551
 552   if(strstr(option,"tim")){
 553     gBenchmark->Stop("EMCALDigitizer");
 554     AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event",
 555          gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
 556   }
 557 }
 558
 559 //____________________________________________________________________________
 560 Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks)
 561 {
 562   //  Returns the shortest time among all time ticks
 563
 564   ticks->Sort() ; //Sort in accordance with times of ticks
 565   TIter it(ticks) ;
 566   AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
 567   Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 568
 569   AliEMCALTick * t ;
 570   while((t=(AliEMCALTick*) it.Next())){
 571     if(t->GetTime() < time)  //This tick starts before crossing
 572       *ctick+=*t ;
 573     else
 574       return time ;
 575
 576     time = ctick->CrossingTime(fTimeThreshold) ;
 577   }
 578   return time ;
 579 }
 580
 581 //____________________________________________________________________________
 582 Bool_t AliEMCALDigitizer::Init()
 583 {
 584   // Makes all memory allocations
 585   fInit = kTRUE ;
 586   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 587
 588   if ( emcalLoader == 0 ) {
 589     Fatal("Init", "Could not obtain the AliEMCALLoader");
 590     return kFALSE;
 591   }
 592
 593   fFirstEvent = 0 ;
 594   fLastEvent = fFirstEvent ;
 595
 596   if(fManager)
 597     fInput = fManager->GetNinputs() ;
 598   else
 599     fInput           = 1 ;
 600
 601   fInputFileNames  = new TString[fInput] ;
 602   fEventNames      = new TString[fInput] ;
 603   fInputFileNames[0] = GetTitle() ;
 604   fEventNames[0]     = fEventFolderName.Data() ;
 605   Int_t index ;
 606   for (index = 1 ; index < fInput ; index++) {
 607     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 608     TString tempo = fManager->GetInputFolderName(index) ;
 609     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager
 610   }
 611
 612   //to prevent cleaning of this object while GetEvent is called
 613   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 614
 615   //PH  Print();
 616   //Calibration instance
 617   fCalibData = emcalLoader->CalibData();
 618   return fInit ;
 619 }
 620
 621 //____________________________________________________________________________
 622 void AliEMCALDigitizer::InitParameters()
 623 {
 624   // Parameter initialization for digitizer
 625   // Tune parameters - 24-nov-04; Apr 29, 2007
 626   // New parameters JLK 14-Apr-2008
 627
 628   fMeanPhotonElectron = 4400;  // electrons per GeV
 629   fPinNoise           = 0.037; // pin noise in GEV from analysis test beam data
 630   if (fPinNoise == 0. )
 631     Warning("InitParameters", "No noise added\n") ;
 632   fDigitThreshold     = fPinNoise * 3; // 3 * sigma
 633   fTimeResolution     = 0.3e-9 ; // 300 psc
 634   fTimeSignalLength   = 1.0e-9 ;
 635
 636   // These defaults are normally not used.
 637   // Values are read from calibration database instead
 638   fADCchannelEC    = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
 639   fADCpedestalEC   = 0.0 ;  // GeV
 640   fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
 641
 642   fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
 643   // hists. for control; no hists on default
 644   fControlHists = 0;
 645   fHists        = 0;
 646 }
 647
 648 //__________________________________________________________________
 649 void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
 650 {
 651   // Allows to produce digits by superimposing background and signal event.
 652   // It is assumed, that headers file with SIGNAL events is opened in
 653   // the constructor.
 654   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed
 655   // Thus we avoid writing (changing) huge and expensive
 656   // backgound files: all output will be writen into SIGNAL, i.e.
 657   // opened in constructor file.
 658   //
 659   // One can open as many files to mix with as one needs.
 660   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
 661   // can be mixed.
 662
 663   if( strcmp(GetName(), "") == 0 )
 664     Init() ;
 665
 666   if(fManager){
 667     Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
 668     return ;
 669   }
 670   // looking for file which contains AliRun
 671   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
 672     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
 673     return ;
 674   }
 675   // looking for the file which contains SDigits
 676   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 677   TString fileName( emcalLoader->GetSDigitsFileName() ) ;
 678     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 679       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
 680     if ( (gSystem->AccessPathName(fileName)) ) {
 681       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
 682       return ;
 683     }
 684     // need to increase the arrays
 685     // MvL: This code only works when fInput == 1, I think.
 686     TString tempo = fInputFileNames[fInput-1] ;
 687     delete [] fInputFileNames ;
 688     fInputFileNames = new TString[fInput+1] ;
 689     fInputFileNames[fInput-1] = tempo ;
 690
 691     tempo = fEventNames[fInput-1] ;
 692     delete [] fEventNames ;
 693     fEventNames = new TString[fInput+1] ;
 694     fEventNames[fInput-1] = tempo ;
 695
 696     fInputFileNames[fInput] = alirunFileName ;
 697     fEventNames[fInput]     = eventFolderName ;
 698     fInput++ ;
 699 }
 700
 701 //__________________________________________________________________
 702 void AliEMCALDigitizer::Print1(Option_t * option)
 703 { // 19-nov-04 - just for convinience
 704   Print();
 705   PrintDigits(option);
 706 }
 707
 708 //__________________________________________________________________
 709 void AliEMCALDigitizer::Print(Option_t*)const
 710 {
 711   // Print Digitizer's parameters
 712   printf("Print: \n------------------- %s -------------", GetName() ) ;
 713   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 714     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 715
 716     Int_t nStreams ;
 717     if (fManager)
 718       nStreams =  GetNInputStreams() ;
 719     else
 720       nStreams = fInput ;
 721
 722     Int_t index = 0 ;
 723
 724     AliRunLoader *rl=0;
 725
 726     for (index = 0 ; index < nStreams ; index++) {
 727       TString tempo(fEventNames[index]) ;
 728       tempo += index ;
 729       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
 730         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
 731       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 732       TString fileName( emcalLoader->GetSDigitsFileName() ) ;
 733       if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 734         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
 735       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
 736     }
 737
 738     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 739
 740     printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
 741
 742     printf("\nWith following parameters:\n") ;
 743
 744     printf("    Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise) ;
 745     printf("    Threshold  in EMC  (fDigitThreshold) = %f\n", fDigitThreshold)  ;
 746     printf("---------------------------------------------------\n")  ;
 747   }
 748   else
 749     printf("Print: AliEMCALDigitizer not initialized") ;
 750 }
 751
 752 //__________________________________________________________________
 753 void AliEMCALDigitizer::PrintDigits(Option_t * option)
 754 {
 755   //utility method for printing digit information
 756
 757   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 758   TClonesArray * digits  = emcalLoader->Digits() ;
 759   TClonesArray * sdigits = emcalLoader->SDigits() ;
 760
 761   printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
 762   printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
 763
 764   if(strstr(option,"all")){
 765     //loop over digits
 766     AliEMCALDigit * digit;
 767     printf("\nEMC digits (with primaries):\n")  ;
 768     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 769     Int_t index ;
 770     for (index = 0 ; index < digits->GetEntries() ; index++) {
 771       digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
 772       printf("\n%6d  %8d    %6.5e %4d      %2d : ",
 773               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 774       Int_t iprimary;
 775       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 776         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 777       }
 778     }
 779   }
 780   printf("\n");
 781 }
 782
 783 //__________________________________________________________________
 784 Float_t AliEMCALDigitizer::TimeOfNoise(void)
 785 {
 786   // Calculates the time signal generated by noise
 787   //PH  Info("TimeOfNoise", "Change me") ;
 788   return gRandom->Rndm() * 1.28E-5;
 789 }
 790
 791 //__________________________________________________________________
 792 void AliEMCALDigitizer::Unload()
 793 {
 794   // Unloads the SDigits and Digits
 795   AliRunLoader *rl=0;
 796
 797   Int_t i ;
 798   for(i = 1 ; i < fInput ; i++){
 799     TString tempo(fEventNames[i]) ;
 800     tempo += i;
 801     if ((rl = AliRunLoader::GetRunLoader(tempo)))
 802       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
 803   }
 804   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 805   emcalLoader->UnloadDigits() ;
 806 }
 807
 808 //_________________________________________________________________________________________
 809 void AliEMCALDigitizer::WriteDigits()
 810 {
 811
 812   // Makes TreeD in the output file.
 813   // Check if branch already exists:
 814   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
 815   //      already existing branches.
 816   //   else creates branch with Digits, named "EMCAL", title "...",
 817   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
 818   //      and names of files, from which digits are made.
 819
 820   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 821
 822   const TClonesArray * digits = emcalLoader->Digits() ;
 823   TTree * treeD = emcalLoader->TreeD();
 824   if ( !treeD ) {
 825     emcalLoader->MakeDigitsContainer();
 826     treeD = emcalLoader->TreeD();
 827   }
 828
 829   // -- create Digits branch
 830   Int_t bufferSize = 32000 ;
 831   TBranch * digitsBranch = 0;
 832   if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
 833     digitsBranch->SetAddress(&digits);
 834     AliWarning("Digits Branch already exists. Not all digits will be visible");
 835   }
 836   else
 837     treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
 838   //digitsBranch->SetTitle(fEventFolderName);
 839   treeD->Fill() ;
 840
 841   emcalLoader->WriteDigits("OVERWRITE");
 842   emcalLoader->WriteDigitizer("OVERWRITE");
 843
 844   Unload() ;
 845
 846 }
 847
 848 //__________________________________________________________________
 849 void AliEMCALDigitizer::Browse(TBrowser* b)
 850 {
 851   if(fHists) b->Add(fHists);
 852   TTask::Browse(b);
 853 }
 854
 855 //__________________________________________________________________
 856 TList *AliEMCALDigitizer::BookControlHists(int var)
 857 {
 858   // 22-nov-04
 859   // histograms for monitoring digitizer performance
 860
 861   Info("BookControlHists"," started ");
 862   gROOT->cd();
 863   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 864   if(var>=1){
 865     new TH1F("hDigiN",  "#EMCAL digits with fAmp > fDigitThreshold",
 866     fNADCEC+1, -0.5, Double_t(fNADCEC));
 867     new TH1F("HDigiSumEnergy","Sum.EMCAL energy from digi", 1000, 0.0, 200.);
 868     new TH1F("hDigiAmp",  "EMCAL digital amplitude", fNADCEC+1, -0.5, Double_t(fNADCEC));
 869     new TH1F("hDigiEnergy","EMCAL cell energy", 2000, 0.0, 200.);
 870     new TH1F("hDigiAbsId","EMCAL absId cells with fAmp > fDigitThreshold ",
 871     geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
 872   }
 873
 874   fHists = AliEMCALHistoUtilities::MoveHistsToList("EmcalDigiControlHists", kFALSE);
 875   fHists = 0; //huh? JLK 03-Mar-2006
 876   return fHists;
 877 }
 878
 879 //__________________________________________________________________
 880 void AliEMCALDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
 881 {
 882   AliEMCALHistoUtilities::SaveListOfHists(fHists, name, kSingleKey, opt);
 883 }