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