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