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)
 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)
 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 isTrd1Geom = -1; // -1 - mean undefined
 245   static int nEMC=0; //max number of digits possible
 246
 247   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 248   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 249   Int_t readEvent = event ;
 250   // fManager is data member from AliDigitizer
 251   if (fManager)
 252     readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
 253   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 254                   readEvent, GetTitle(), fEventFolderName.Data())) ;
 255   rl->GetEvent(readEvent);
 256
 257   TClonesArray * digits = emcalLoader->Digits() ;
 258   digits->Delete() ;
 259
 260   // Load Geometry
 261   AliEMCALGeometry *geom = 0;
 262   if (rl->GetAliRun()) {
 263     AliEMCAL * emcal  = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
 264     geom = emcal->GetGeometry();
 265   }
 266   else
 267     AliFatal("Could not get AliRun from runLoader");
 268
 269   if(isTrd1Geom < 0) {
 270     AliInfo(Form(" get Geometry %s : %s ", geom->GetName(),geom->GetTitle()));
 271     TString ng(geom->GetName());
 272     isTrd1Geom = 0;
 273     if(ng.Contains("SHISH") &&  ng.Contains("TRD1")) isTrd1Geom = 1;
 274
 275     if(isTrd1Geom == 0) nEMC = geom->GetNPhi()*geom->GetNZ();
 276     else                nEMC = geom->GetNCells();
 277     AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s | isTrd1Geom %i\n", nEMC, geom->GetName(), isTrd1Geom));
 278   }
 279   Int_t absID ;
 280
 281   digits->Expand(nEMC) ;
 282
 283   // get first the sdigitizer from the tasks list (must have same name as the digitizer)
 284   if ( !emcalLoader->SDigitizer() )
 285     emcalLoader->LoadSDigitizer();
 286   AliEMCALSDigitizer * sDigitizer = dynamic_cast<AliEMCALSDigitizer*>(emcalLoader->SDigitizer());
 287
 288   if ( !sDigitizer )
 289     Fatal("Digitize", "SDigitizer with name %s %s not found", fEventFolderName.Data(), GetTitle() ) ;
 290
 291   //take all the inputs to add together and load the SDigits
 292   TObjArray * sdigArray = new TObjArray(fInput) ;
 293   sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
 294   Int_t i ;
 295
 296   for(i = 1 ; i < fInput ; i++){
 297     TString tempo(fEventNames[i]) ;
 298     tempo += i ;
 299
 300     AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ;
 301
 302     if (rl2==0)
 303       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
 304
 305     if (fManager)
 306       readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
 307     Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
 308     rl2->LoadSDigits();
 309     rl2->GetEvent(readEvent);
 310     AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
 311     sdigArray->AddAt(emcalLoader2->SDigits(), i) ;
 312   }
 313
 314   //Find the first tower with signal
 315   Int_t nextSig = nEMC + 1 ;
 316   TClonesArray * sdigits ;
 317   for(i = 0 ; i < fInput ; i++){
 318     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 319     if ( !sdigits->GetEntriesFast() )
 320       continue ;
 321     Int_t curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(0))->GetId() ;
 322      if(curNext < nextSig)
 323        nextSig = curNext ;
 324      AliDebug(1,Form("input %i : #sdigits %i \n",
 325                      i, sdigits->GetEntriesFast()));
 326   }
 327   AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
 328
 329   TArrayI index(fInput) ;
 330   index.Reset() ;  //Set all indexes to zero
 331
 332   AliEMCALDigit * digit ;
 333   AliEMCALDigit * curSDigit ;
 334
 335   TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ;
 336
 337   //Put Noise contribution
 338   for(absID = 0; absID < nEMC; absID++){ // Nov 30, 2006 by PAI; was from 1 to nEMC
 339     Float_t amp = 0 ;
 340     // amplitude set to zero, noise will be added later
 341     new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0, TimeOfNoise() ); // absID-1->absID
 342     //look if we have to add signal?
 343     digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
 344
 345     if(absID==nextSig){
 346       //Add SDigits from all inputs
 347       ticks->Clear() ;
 348       Int_t contrib = 0 ;
 349       Float_t a = digit->GetAmp() ;
 350       Float_t b = TMath::Abs( a /fTimeSignalLength) ;
 351       //Mark the beginning of the signal
 352       new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime(),0, b);
 353       //Mark the end of the signal
 354       new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime()+fTimeSignalLength, -a, -b);
 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           a = curSDigit->GetAmp() ;
 373           b = a /fTimeSignalLength ;
 374           new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime(),0, b);
 375           new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
 376
 377           *digit = *digit + *curSDigit ;  //add energies
 378
 379           index[i]++ ;
 380           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 381             curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 382           else
 383             curSDigit = 0 ;
 384         }
 385       }
 386       // add fluctuations for photo-electron creation
 387       amp = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
 388       amp *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
 389
 390       //calculate and set time
 391       Float_t time = FrontEdgeTime(ticks) ;
 392       digit->SetTime(time) ;
 393
 394       //Find next signal module
 395       nextSig = nEMC + 1 ;
 396       for(i = 0 ; i < fInput ; i++){
 397         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 398         Int_t curNext = nextSig ;
 399         if(sdigits->GetEntriesFast() > index[i] ){
 400           curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]))->GetId() ;
 401         }
 402         if(curNext < nextSig) nextSig = curNext ;
 403       }
 404     }
 405     // add the noise now
 406     amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
 407     digit->SetAmp(sDigitizer->Digitize(amp)) ;
 408     AliDebug(10,Form(" absID %5i amp %f nextSig %5i\n",
 409                      absID, amp, nextSig));
 410   } // for(absID = 1; absID <= nEMC; absID++)
 411
 412   ticks->Delete() ;
 413   delete ticks ;
 414
 415   delete sdigArray ; //We should not delete its contents
 416
 417   //remove digits below thresholds
 418   for(i = 0 ; i < nEMC ; i++){
 419     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
 420     Float_t threshold = fDigitThreshold ;
 421     if(sDigitizer->Calibrate( digit->GetAmp() ) < threshold)
 422       digits->RemoveAt(i) ;
 423     else
 424       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
 425   }
 426
 427   digits->Compress() ;
 428
 429   Int_t ndigits = digits->GetEntriesFast() ;
 430
 431   //Set indexes in list of digits and fill hists.
 432   AliEMCALHistoUtilities::FillH1(fHists, 0, Double_t(ndigits));
 433   Float_t energy=0., esum=0.;
 434   for (i = 0 ; i < ndigits ; i++) {
 435     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
 436     digit->SetIndexInList(i) ;
 437     energy = sDigitizer->Calibrate(digit->GetAmp()) ;
 438     esum += energy;
 439     digit->SetAmp(DigitizeEnergy(energy, digit->GetId()) ) ; // for what ??
 440     AliEMCALHistoUtilities::FillH1(fHists, 2, double(digit->GetAmp()));
 441     AliEMCALHistoUtilities::FillH1(fHists, 3, double(energy));
 442     AliEMCALHistoUtilities::FillH1(fHists, 4, double(digit->GetId()));
 443     //    if(digit->GetId() == nEMC) {
 444     //  printf(" i %i \n", i );
 445     //  digit->Dump();
 446     //  assert(0);
 447     //}
 448   }
 449   AliEMCALHistoUtilities::FillH1(fHists, 1, esum);
 450 }
 451
 452 // //_____________________________________________________________________
 453 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
 454 {
 455   // Returns digitized value of the energy in a cell absId
 456
 457   // Load Geometry
 458   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 459
 460   if (geom==0)
 461     AliFatal("Did not get geometry from EMCALLoader");
 462
 463   Int_t iSupMod = -1;
 464   Int_t nModule  = -1;
 465   Int_t nIphi   = -1;
 466   Int_t nIeta   = -1;
 467   Int_t iphi    = -1;
 468   Int_t ieta    = -1;
 469   Int_t channel = -999;
 470
 471   Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
 472   if(!bCell)
 473     Error("DigitizeEnergy","Wrong cell id number : AbsId %i ", AbsId) ;
 474   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 475
 476   if(fCalibData) {
 477     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 478     fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
 479   }
 480
 481   channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
 482
 483   if(channel > fNADCEC )
 484     channel =  fNADCEC ;
 485   return channel ;
 486
 487 }
 488
 489 //____________________________________________________________________________
 490 void AliEMCALDigitizer::Exec(Option_t *option)
 491 {
 492   // Steering method to process digitization for events
 493   // in the range from fFirstEvent to fLastEvent.
 494   // This range is optionally set by SetEventRange().
 495   // if fLastEvent=-1, then process events until the end.
 496   // by default fLastEvent = fFirstEvent (process only one event)
 497
 498   if (!fInit) { // to prevent overwrite existing file
 499     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
 500     return ;
 501   }
 502
 503   if (strstr(option,"print")) {
 504
 505     Print();
 506     return ;
 507   }
 508
 509   if(strstr(option,"tim"))
 510     gBenchmark->Start("EMCALDigitizer");
 511
 512   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 513   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 514
 515   // Post Digitizer to the white board
 516   emcalLoader->PostDigitizer(this) ;
 517
 518   if (fLastEvent == -1)  {
 519     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 520   }
 521   else if (fManager)
 522     fLastEvent = fFirstEvent ; // what is this ??
 523
 524   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 525   Int_t ievent;
 526
 527   rl->LoadSDigits("EMCAL");
 528   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 529
 530     rl->GetEvent(ievent);
 531
 532     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 533
 534     WriteDigits() ;
 535
 536     if(strstr(option,"deb"))
 537       PrintDigits(option);
 538     if(strstr(option,"table")) gObjectTable->Print();
 539
 540     //increment the total number of Digits per run
 541     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
 542   }
 543
 544   emcalLoader->CleanDigitizer() ;
 545
 546   if(strstr(option,"tim")){
 547     gBenchmark->Stop("EMCALDigitizer");
 548     AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event",
 549          gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
 550   }
 551 }
 552
 553 //____________________________________________________________________________
 554 Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks)
 555 {
 556   //  Returns the shortest time among all time ticks
 557
 558   ticks->Sort() ; //Sort in accordance with times of ticks
 559   TIter it(ticks) ;
 560   AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
 561   Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 562
 563   AliEMCALTick * t ;
 564   while((t=(AliEMCALTick*) it.Next())){
 565     if(t->GetTime() < time)  //This tick starts before crossing
 566       *ctick+=*t ;
 567     else
 568       return time ;
 569
 570     time = ctick->CrossingTime(fTimeThreshold) ;
 571   }
 572   return time ;
 573 }
 574
 575 //____________________________________________________________________________
 576 Bool_t AliEMCALDigitizer::Init()
 577 {
 578   // Makes all memory allocations
 579   fInit = kTRUE ;
 580   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 581
 582   if ( emcalLoader == 0 ) {
 583     Fatal("Init", "Could not obtain the AliEMCALLoader");
 584     return kFALSE;
 585   }
 586
 587   fFirstEvent = 0 ;
 588   fLastEvent = fFirstEvent ;
 589
 590   if(fManager)
 591     fInput = fManager->GetNinputs() ;
 592   else
 593     fInput           = 1 ;
 594
 595   fInputFileNames  = new TString[fInput] ;
 596   fEventNames      = new TString[fInput] ;
 597   fInputFileNames[0] = GetTitle() ;
 598   fEventNames[0]     = fEventFolderName.Data() ;
 599   Int_t index ;
 600   for (index = 1 ; index < fInput ; index++) {
 601     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 602     TString tempo = fManager->GetInputFolderName(index) ;
 603     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager
 604   }
 605
 606   //to prevent cleaning of this object while GetEvent is called
 607   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 608
 609   //PH  Print();
 610   //Calibration instance
 611   fCalibData = emcalLoader->CalibData();
 612   return fInit ;
 613 }
 614
 615 //____________________________________________________________________________
 616 void AliEMCALDigitizer::InitParameters()
 617 {
 618   // Parameter initialization for digitizer
 619   // Tune parameters - 24-nov-04; Apr 29, 2007
 620
 621   fMeanPhotonElectron = 3300;  // electrons per GeV
 622   fPinNoise           = 0.010; // pin noise in GEV from analysis test beam data
 623   if (fPinNoise == 0. )
 624     Warning("InitParameters", "No noise added\n") ;
 625   fDigitThreshold     = fPinNoise * 3; // 3 * sigma
 626   fTimeResolution     = 0.3e-9 ; // 300 psc
 627   fTimeSignalLength   = 1.0e-9 ;
 628
 629   // These defaults are normally not used.
 630   // Values are read from calibration database instead
 631   fADCchannelEC    = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
 632   fADCpedestalEC   = 0.0 ;  // GeV
 633   fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
 634
 635   fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
 636   // hists. for control; no hists on default
 637   fControlHists = 0;
 638   fHists        = 0;
 639 }
 640
 641 //__________________________________________________________________
 642 void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
 643 {
 644   // Allows to produce digits by superimposing background and signal event.
 645   // It is assumed, that headers file with SIGNAL events is opened in
 646   // the constructor.
 647   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed
 648   // Thus we avoid writing (changing) huge and expensive
 649   // backgound files: all output will be writen into SIGNAL, i.e.
 650   // opened in constructor file.
 651   //
 652   // One can open as many files to mix with as one needs.
 653   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
 654   // can be mixed.
 655
 656   if( strcmp(GetName(), "") == 0 )
 657     Init() ;
 658
 659   if(fManager){
 660     Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
 661     return ;
 662   }
 663   // looking for file which contains AliRun
 664   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
 665     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
 666     return ;
 667   }
 668   // looking for the file which contains SDigits
 669   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 670   TString fileName( emcalLoader->GetSDigitsFileName() ) ;
 671     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 672       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
 673     if ( (gSystem->AccessPathName(fileName)) ) {
 674       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
 675       return ;
 676     }
 677     // need to increase the arrays
 678     // MvL: This code only works when fInput == 1, I think.
 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     AliWarning("Digits Branch already exists. Not all digits will be visible");
 827   }
 828   else
 829     treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
 830   //digitsBranch->SetTitle(fEventFolderName);
 831   treeD->Fill() ;
 832
 833   emcalLoader->WriteDigits("OVERWRITE");
 834   emcalLoader->WriteDigitizer("OVERWRITE");
 835
 836   Unload() ;
 837
 838 }
 839
 840 void AliEMCALDigitizer::Browse(TBrowser* b)
 841 {
 842   if(fHists) b->Add(fHists);
 843   TTask::Browse(b);
 844 }
 845
 846 TList *AliEMCALDigitizer::BookControlHists(int var)
 847 {
 848   // 22-nov-04
 849   // histograms for monitoring digitizer performance
 850
 851   Info("BookControlHists"," started ");
 852   gROOT->cd();
 853   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 854   if(var>=1){
 855     new TH1F("hDigiN",  "#EMCAL digits with fAmp > fDigitThreshold",
 856     fNADCEC+1, -0.5, Double_t(fNADCEC));
 857     new TH1F("HDigiSumEnergy","Sum.EMCAL energy from digi", 1000, 0.0, 200.);
 858     new TH1F("hDigiAmp",  "EMCAL digital amplitude", fNADCEC+1, -0.5, Double_t(fNADCEC));
 859     new TH1F("hDigiEnergy","EMCAL cell energy", 2000, 0.0, 200.);
 860     new TH1F("hDigiAbsId","EMCAL absId cells with fAmp > fDigitThreshold ",
 861     geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
 862   }
 863
 864   fHists = AliEMCALHistoUtilities::MoveHistsToList("EmcalDigiControlHists", kFALSE);
 865   fHists = 0; //huh? JLK 03-Mar-2006
 866   return fHists;
 867 }
 868
 869 void AliEMCALDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
 870 {
 871   AliEMCALHistoUtilities::SaveListOfHists(fHists, name, kSingleKey, opt);
 872 }