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