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