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