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; Apr 29, 2007
 619
 620   fMeanPhotonElectron = 3300;  // electrons per GeV
 621   fPinNoise           = 0.010; // pin noise in GEV from analysis test beam data
 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   // These defaults are normally not used.
 629   // Values are read from calibration database instead
 630   fADCchannelEC    = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
 631   fADCpedestalEC   = 0.0 ;  // GeV
 632   fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
 633
 634   fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
 635   // hists. for control; no hists on default
 636   fControlHists = 0;
 637   fHists        = 0;
 638 }
 639
 640 //__________________________________________________________________
 641 void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
 642 {
 643   // Allows to produce digits by superimposing background and signal event.
 644   // It is assumed, that headers file with SIGNAL events is opened in
 645   // the constructor.
 646   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed
 647   // Thus we avoid writing (changing) huge and expensive
 648   // backgound files: all output will be writen into SIGNAL, i.e.
 649   // opened in constructor file.
 650   //
 651   // One can open as many files to mix with as one needs.
 652   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
 653   // can be mixed.
 654
 655   if( strcmp(GetName(), "") == 0 )
 656     Init() ;
 657
 658   if(fManager){
 659     Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
 660     return ;
 661   }
 662   // looking for file which contains AliRun
 663   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
 664     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
 665     return ;
 666   }
 667   // looking for the file which contains SDigits
 668   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 669   TString fileName( emcalLoader->GetSDigitsFileName() ) ;
 670     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 671       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
 672     if ( (gSystem->AccessPathName(fileName)) ) {
 673       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
 674       return ;
 675     }
 676     // need to increase the arrays
 677     // MvL: This code only works when fInput == 1, I think.
 678     TString tempo = fInputFileNames[fInput-1] ;
 679     delete [] fInputFileNames ;
 680     fInputFileNames = new TString[fInput+1] ;
 681     fInputFileNames[fInput-1] = tempo ;
 682
 683     tempo = fEventNames[fInput-1] ;
 684     delete [] fEventNames ;
 685     fEventNames = new TString[fInput+1] ;
 686     fEventNames[fInput-1] = tempo ;
 687
 688     fInputFileNames[fInput] = alirunFileName ;
 689     fEventNames[fInput]     = eventFolderName ;
 690     fInput++ ;
 691 }
 692
 693 void AliEMCALDigitizer::Print1(Option_t * option)
 694 { // 19-nov-04 - just for convinience
 695   Print();
 696   PrintDigits(option);
 697 }
 698
 699 //__________________________________________________________________
 700 void AliEMCALDigitizer::Print(Option_t*)const
 701 {
 702   // Print Digitizer's parameters
 703   printf("Print: \n------------------- %s -------------", GetName() ) ;
 704   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 705     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 706
 707     Int_t nStreams ;
 708     if (fManager)
 709       nStreams =  GetNInputStreams() ;
 710     else
 711       nStreams = fInput ;
 712
 713     Int_t index = 0 ;
 714
 715     AliRunLoader *rl=0;
 716
 717     for (index = 0 ; index < nStreams ; index++) {
 718       TString tempo(fEventNames[index]) ;
 719       tempo += index ;
 720       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
 721         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
 722       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 723       TString fileName( emcalLoader->GetSDigitsFileName() ) ;
 724       if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 725         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
 726       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
 727     }
 728
 729     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 730
 731     printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
 732
 733     printf("\nWith following parameters:\n") ;
 734
 735     printf("    Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise) ;
 736     printf("    Threshold  in EMC  (fDigitThreshold) = %f\n", fDigitThreshold)  ;
 737     printf("---------------------------------------------------\n")  ;
 738   }
 739   else
 740     printf("Print: AliEMCALDigitizer not initialized") ;
 741 }
 742
 743 //__________________________________________________________________
 744 void AliEMCALDigitizer::PrintDigits(Option_t * option)
 745 {
 746   //utility method for printing digit information
 747
 748   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 749   TClonesArray * digits  = emcalLoader->Digits() ;
 750   TClonesArray * sdigits = emcalLoader->SDigits() ;
 751
 752   printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
 753   printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
 754
 755   if(strstr(option,"all")){
 756     //loop over digits
 757     AliEMCALDigit * digit;
 758     printf("\nEMC digits (with primaries):\n")  ;
 759     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 760     Int_t index ;
 761     for (index = 0 ; index < digits->GetEntries() ; index++) {
 762       digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
 763       printf("\n%6d  %8d    %6.5e %4d      %2d : ",
 764               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 765       Int_t iprimary;
 766       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 767         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 768       }
 769     }
 770   }
 771   printf("\n");
 772 }
 773
 774 //__________________________________________________________________
 775 Float_t AliEMCALDigitizer::TimeOfNoise(void)
 776 {
 777   // Calculates the time signal generated by noise
 778   //PH  Info("TimeOfNoise", "Change me") ;
 779   return gRandom->Rndm() * 1.28E-5;
 780 }
 781
 782 //__________________________________________________________________
 783 void AliEMCALDigitizer::Unload()
 784 {
 785   // Unloads the SDigits and Digits
 786   AliRunLoader *rl=0;
 787
 788   Int_t i ;
 789   for(i = 1 ; i < fInput ; i++){
 790     TString tempo(fEventNames[i]) ;
 791     tempo += i;
 792     if ((rl = AliRunLoader::GetRunLoader(tempo)))
 793       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
 794   }
 795   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 796   emcalLoader->UnloadDigits() ;
 797 }
 798
 799 //_________________________________________________________________________________________
 800 void AliEMCALDigitizer::WriteDigits()
 801 {
 802
 803   // Makes TreeD in the output file.
 804   // Check if branch already exists:
 805   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
 806   //      already existing branches.
 807   //   else creates branch with Digits, named "EMCAL", title "...",
 808   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
 809   //      and names of files, from which digits are made.
 810
 811   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 812
 813   const TClonesArray * digits = emcalLoader->Digits() ;
 814   TTree * treeD = emcalLoader->TreeD();
 815   if ( !treeD ) {
 816     emcalLoader->MakeDigitsContainer();
 817     treeD = emcalLoader->TreeD();
 818   }
 819
 820   // -- create Digits branch
 821   Int_t bufferSize = 32000 ;
 822   TBranch * digitsBranch = 0;
 823   if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
 824     digitsBranch->SetAddress(&digits);
 825     AliWarning("Digits Branch already exists. Not all digits will be visible");
 826   }
 827   else
 828     treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
 829   //digitsBranch->SetTitle(fEventFolderName);
 830   treeD->Fill() ;
 831
 832   emcalLoader->WriteDigits("OVERWRITE");
 833   emcalLoader->WriteDigitizer("OVERWRITE");
 834
 835   Unload() ;
 836
 837 }
 838
 839 void AliEMCALDigitizer::Browse(TBrowser* b)
 840 {
 841   if(fHists) b->Add(fHists);
 842   TTask::Browse(b);
 843 }
 844
 845 TList *AliEMCALDigitizer::BookControlHists(int var)
 846 {
 847   // 22-nov-04
 848   // histograms for monitoring digitizer performance
 849
 850   Info("BookControlHists"," started ");
 851   gROOT->cd();
 852   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 853   if(var>=1){
 854     new TH1F("hDigiN",  "#EMCAL digits with fAmp > fDigitThreshold",
 855     fNADCEC+1, -0.5, Double_t(fNADCEC));
 856     new TH1F("HDigiSumEnergy","Sum.EMCAL energy from digi", 1000, 0.0, 200.);
 857     new TH1F("hDigiAmp",  "EMCAL digital amplitude", fNADCEC+1, -0.5, Double_t(fNADCEC));
 858     new TH1F("hDigiEnergy","EMCAL cell energy", 2000, 0.0, 200.);
 859     new TH1F("hDigiAbsId","EMCAL absId cells with fAmp > fDigitThreshold ",
 860     geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
 861   }
 862
 863   fHists = AliEMCALHistoUtilities::MoveHistsToList("EmcalDigiControlHists", kFALSE);
 864   fHists = 0; //huh? JLK 03-Mar-2006
 865   return fHists;
 866 }
 867
 868 void AliEMCALDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
 869 {
 870   AliEMCALHistoUtilities::SaveListOfHists(fHists, name, kSingleKey, opt);
 871 }