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 "AliEMCALDigit.h"
  78 #include "AliEMCAL.h"
  79 #include "AliEMCALLoader.h"
  80 #include "AliEMCALDigitizer.h"
  81 #include "AliEMCALSDigitizer.h"
  82 #include "AliEMCALGeometry.h"
  83 #include "AliEMCALTick.h"
  84 #include "AliEMCALHistoUtilities.h"
  85
  86 ClassImp(AliEMCALDigitizer)
  87
  88
  89 //____________________________________________________________________________
  90 AliEMCALDigitizer::AliEMCALDigitizer()
  91   : AliDigitizer("",""),
  92     fDefaultInit(kTRUE),
  93     fDigitsInRun(0),
  94     fInit(0),
  95     fInput(0),
  96     fInputFileNames(0x0),
  97     fEventNames(0x0),
  98     fDigitThreshold(0),
  99     fMeanPhotonElectron(0),
 100     fPedestal(0),
 101     fSlope(0),
 102     fPinNoise(0),
 103     fTimeResolution(0),
 104     fTimeThreshold(0),
 105     fTimeSignalLength(0),
 106     fADCchannelEC(0),
 107     fADCpedestalEC(0),
 108     fNADCEC(0),
 109     fEventFolderName(""),
 110     fFirstEvent(0),
 111     fLastEvent(0),
 112     fControlHists(0),
 113     fHists(0)
 114 {
 115   // ctor
 116   InitParameters() ;
 117   fManager = 0 ;                     // We work in the standalong mode
 118 }
 119
 120 //____________________________________________________________________________
 121 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
 122   : AliDigitizer("EMCAL"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
 123     fDefaultInit(kFALSE),
 124     fDigitsInRun(0),
 125     fInit(0),
 126     fInput(0),
 127     fInputFileNames(0),
 128     fEventNames(0),
 129     fDigitThreshold(0),
 130     fMeanPhotonElectron(0),
 131     fPedestal(0),
 132     fSlope(0),
 133     fPinNoise(0),
 134     fTimeResolution(0),
 135     fTimeThreshold(0),
 136     fTimeSignalLength(0),
 137     fADCchannelEC(0),
 138     fADCpedestalEC(0),
 139     fNADCEC(0),
 140     fEventFolderName(eventFolderName),
 141     fFirstEvent(0),
 142     fLastEvent(0),
 143     fControlHists(0),
 144     fHists(0)
 145 {
 146   // ctor
 147   InitParameters() ;
 148   Init() ;
 149   fManager = 0 ;                     // We work in the standalong mode
 150 }
 151
 152 //____________________________________________________________________________
 153 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
 154   : AliDigitizer(d.GetName(),d.GetTitle()),
 155     fDefaultInit(d.fDefaultInit),
 156     fDigitsInRun(d.fDigitsInRun),
 157     fInit(d.fInit),
 158     fInput(d.fInput),
 159     fInputFileNames(d.fInputFileNames),
 160     fEventNames(d.fEventNames),
 161     fDigitThreshold(d.fDigitThreshold),
 162     fMeanPhotonElectron(d.fMeanPhotonElectron),
 163     fPedestal(d.fPedestal),
 164     fSlope(d.fSlope),
 165     fPinNoise(d.fPinNoise),
 166     fTimeResolution(d.fTimeResolution),
 167     fTimeThreshold(d.fTimeThreshold),
 168     fTimeSignalLength(d.fTimeSignalLength),
 169     fADCchannelEC(d.fADCchannelEC),
 170     fADCpedestalEC(d.fADCpedestalEC),
 171     fNADCEC(d.fNADCEC),
 172     fEventFolderName(d.fEventFolderName),
 173     fFirstEvent(d.fFirstEvent),
 174     fLastEvent(d.fLastEvent),
 175     fControlHists(d.fControlHists),
 176     fHists(d.fHists)
 177 {
 178   // copyy ctor
 179  }
 180
 181 //____________________________________________________________________________
 182 AliEMCALDigitizer::AliEMCALDigitizer(AliRunDigitizer * rd)
 183   : AliDigitizer(rd,"EMCAL"+AliConfig::Instance()->GetDigitizerTaskName()),
 184     fDefaultInit(kFALSE),
 185     fDigitsInRun(0),
 186     fInit(0),
 187     fInput(0),
 188     fInputFileNames(0),
 189     fEventNames(0),
 190     fDigitThreshold(0.),
 191     fMeanPhotonElectron(0),
 192     fPedestal(0),
 193     fSlope(0.),
 194     fPinNoise(0),
 195     fTimeResolution(0.),
 196     fTimeThreshold(0),
 197     fTimeSignalLength(0),
 198     fADCchannelEC(0),
 199     fADCpedestalEC(0),
 200     fNADCEC(0),
 201     fEventFolderName(0),
 202     fFirstEvent(0),
 203     fLastEvent(0),
 204     fControlHists(0),
 205     fHists(0)
 206 {
 207   // ctor Init() is called by RunDigitizer
 208   fManager = rd ;
 209   fEventFolderName = fManager->GetInputFolderName(0) ;
 210   SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
 211   InitParameters() ;
 212 }
 213
 214 //____________________________________________________________________________
 215   AliEMCALDigitizer::~AliEMCALDigitizer()
 216 {
 217   //dtor
 218   if (AliRunLoader::GetRunLoader()) {
 219     AliLoader *emcalLoader=0;
 220     if ((emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")))
 221       emcalLoader->CleanDigitizer();
 222   }
 223   else
 224     AliDebug(1," no runloader present");
 225   delete [] fInputFileNames ;
 226   delete [] fEventNames ;
 227
 228   if(fHists) delete fHists;
 229 }
 230
 231 //____________________________________________________________________________
 232 void AliEMCALDigitizer::Digitize(Int_t event)
 233 {
 234
 235   // Makes the digitization of the collected summable digits
 236   // for this it first creates the array of all EMCAL modules
 237   // filled with noise (different for EMC, CPV and PPSD) and
 238   // after that adds contributions from SDigits. This design
 239   // helps to avoid scanning over the list of digits to add
 240   // contribution of any new SDigit.
 241   static int isTrd1Geom = -1; // -1 - mean undefined
 242   static int nEMC=0; //max number of digits possible
 243
 244   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 245   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 246   Int_t readEvent = event ;
 247   // fManager is data member from AliDigitizer
 248   if (fManager)
 249     readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
 250   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 251                   readEvent, GetTitle(), fEventFolderName.Data())) ;
 252   rl->GetEvent(readEvent);
 253
 254   TClonesArray * digits = emcalLoader->Digits() ;
 255   digits->Clear() ;
 256
 257   // Load Geometry
 258   // const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 259   rl->LoadgAlice();
 260   AliRun * gAlice = rl->GetAliRun();
 261   AliEMCAL * emcal  = (AliEMCAL*)gAlice->GetDetector("EMCAL");
 262   AliEMCALGeometry * geom = emcal->GetGeometry();
 263
 264   if(isTrd1Geom < 0) {
 265     TString ng(geom->GetName());
 266     isTrd1Geom = 0;
 267     if(ng.Contains("SHISH") &&  ng.Contains("TRD1")) isTrd1Geom = 1;
 268
 269     if(isTrd1Geom == 0) nEMC = geom->GetNPhi()*geom->GetNZ();
 270     else                nEMC = geom->GetNCells();
 271     AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s | isTrd1Geom %i\n", nEMC, geom->GetName(), isTrd1Geom));
 272   }
 273   Int_t absID ;
 274
 275   digits->Expand(nEMC) ;
 276
 277   // get first the sdigitizer from the tasks list (must have same name as the digitizer)
 278   if ( !emcalLoader->SDigitizer() )
 279     emcalLoader->LoadSDigitizer();
 280   AliEMCALSDigitizer * sDigitizer = dynamic_cast<AliEMCALSDigitizer*>(emcalLoader->SDigitizer());
 281
 282   if ( !sDigitizer )
 283     Fatal("Digitize", "SDigitizer with name %s %s not found", fEventFolderName.Data(), GetTitle() ) ;
 284
 285   //take all the inputs to add together and load the SDigits
 286   TObjArray * sdigArray = new TObjArray(fInput) ;
 287   sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
 288   Int_t i ;
 289
 290   for(i = 1 ; i < fInput ; i++){
 291     TString tempo(fEventNames[i]) ;
 292     tempo += i ;
 293
 294     AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ;
 295
 296     if (rl2==0)
 297       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
 298
 299     if (fManager)
 300       readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
 301     Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
 302     rl2->LoadSDigits();
 303     rl2->GetEvent(readEvent);
 304     AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
 305     sdigArray->AddAt(emcalLoader2->SDigits(), i) ;
 306   }
 307
 308   //Find the first tower with signal
 309   Int_t nextSig = nEMC + 1 ;
 310   TClonesArray * sdigits ;
 311   for(i = 0 ; i < fInput ; i++){
 312     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 313     if ( !sdigits->GetEntriesFast() )
 314       continue ;
 315     Int_t curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(0))->GetId() ;
 316      if(curNext < nextSig)
 317        nextSig = curNext ;
 318      AliDebug(1,Form("input %i : #sdigits %i \n",
 319                      i, sdigits->GetEntriesFast()));
 320   }
 321   AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
 322
 323   TArrayI index(fInput) ;
 324   index.Reset() ;  //Set all indexes to zero
 325
 326   AliEMCALDigit * digit ;
 327   AliEMCALDigit * curSDigit ;
 328
 329   TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ;
 330
 331   //Put Noise contribution
 332   for(absID = 1; absID <= nEMC; absID++){
 333     Float_t amp = 0 ;
 334     // amplitude set to zero, noise will be added later
 335     new((*digits)[absID-1]) AliEMCALDigit( -1, -1, absID, 0, TimeOfNoise() ) ;
 336     //look if we have to add signal?
 337     digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID-1)) ;
 338
 339     if(absID==nextSig){
 340       //Add SDigits from all inputs
 341       ticks->Clear() ;
 342       Int_t contrib = 0 ;
 343       Float_t a = digit->GetAmp() ;
 344       Float_t b = TMath::Abs( a /fTimeSignalLength) ;
 345       //Mark the beginning of the signal
 346       new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime(),0, b);
 347       //Mark the end of the signal
 348       new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime()+fTimeSignalLength, -a, -b);
 349
 350       // loop over input
 351       for(i = 0; i< fInput ; i++){  //loop over (possible) merge sources
 352         if(dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 353           curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 354         else
 355           curSDigit = 0 ;
 356         //May be several digits will contribute from the same input
 357         while(curSDigit && (curSDigit->GetId() == absID)){
 358           //Shift primary to separate primaries belonging different inputs
 359           Int_t primaryoffset ;
 360           if(fManager)
 361             primaryoffset = fManager->GetMask(i) ;
 362           else
 363             primaryoffset = i ;
 364           curSDigit->ShiftPrimary(primaryoffset) ;
 365
 366           a = curSDigit->GetAmp() ;
 367           b = a /fTimeSignalLength ;
 368           new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime(),0, b);
 369           new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
 370
 371           *digit = *digit + *curSDigit ;  //add energies
 372
 373           index[i]++ ;
 374           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 375             curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 376           else
 377             curSDigit = 0 ;
 378         }
 379       }
 380       // add fluctuations for photo-electron creation
 381       amp = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
 382       amp *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
 383
 384       //calculate and set time
 385       Float_t time = FrontEdgeTime(ticks) ;
 386       digit->SetTime(time) ;
 387
 388       //Find next signal module
 389       nextSig = nEMC + 1 ;
 390       for(i = 0 ; i < fInput ; i++){
 391         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 392         Int_t curNext = nextSig ;
 393         if(sdigits->GetEntriesFast() > index[i] ){
 394           curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]))->GetId() ;
 395         }
 396         if(curNext < nextSig) nextSig = curNext ;
 397       }
 398     }
 399     // add the noise now
 400     amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
 401     digit->SetAmp(sDigitizer->Digitize(amp)) ;
 402     AliDebug(10,Form(" absID %5i amp %f nextSig %5i\n",
 403                      absID, amp, nextSig));
 404   } // for(absID = 1; absID <= nEMC; absID++)
 405
 406   ticks->Delete() ;
 407   delete ticks ;
 408
 409   delete sdigArray ; //We should not delete its contents
 410
 411   //remove digits below thresholds
 412   for(i = 0 ; i < nEMC ; i++){
 413     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
 414     Float_t threshold = fDigitThreshold ;
 415     if(sDigitizer->Calibrate( digit->GetAmp() ) < threshold)
 416       digits->RemoveAt(i) ;
 417     else
 418       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
 419   }
 420
 421   digits->Compress() ;
 422
 423   Int_t ndigits = digits->GetEntriesFast() ;
 424   digits->Expand(ndigits) ;
 425
 426   //Set indexes in list of digits and fill hists.
 427   AliEMCALHistoUtilities::FillH1(fHists, 0, Double_t(ndigits));
 428   Float_t energy=0., esum=0.;
 429   for (i = 0 ; i < ndigits ; i++) {
 430     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
 431     digit->SetIndexInList(i) ;
 432     energy = sDigitizer->Calibrate(digit->GetAmp()) ;
 433     esum += energy;
 434     digit->SetAmp(DigitizeEnergy(energy, digit->GetId()) ) ; // for what ??
 435     AliEMCALHistoUtilities::FillH1(fHists, 2, double(digit->GetAmp()));
 436     AliEMCALHistoUtilities::FillH1(fHists, 3, double(energy));
 437     AliEMCALHistoUtilities::FillH1(fHists, 4, double(digit->GetId()));
 438   }
 439   AliEMCALHistoUtilities::FillH1(fHists, 1, esum);
 440 }
 441
 442 // //_____________________________________________________________________
 443 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
 444 {
 445   // Returns digitized value of the energy in a cell absId
 446   // Loader
 447   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 448   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
 449     (rl->GetDetectorLoader("EMCAL"));
 450
 451   // Load EMCAL Geometry
 452   rl->LoadgAlice();
 453   AliRun * gAlice = rl->GetAliRun();
 454   AliEMCAL * emcal  = (AliEMCAL*)gAlice->GetDetector("EMCAL");
 455   AliEMCALGeometry * geom = emcal->GetGeometry();
 456
 457   if (geom==0)
 458     AliFatal("Did not get geometry from EMCALLoader") ;
 459
 460   Int_t iSupMod = -1;
 461   Int_t nTower  = -1;
 462   Int_t nIphi   = -1;
 463   Int_t nIeta   = -1;
 464   Int_t iphi    = -1;
 465   Int_t ieta    = -1;
 466   Int_t channel = -999;
 467
 468   Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nTower, nIphi, nIeta) ;
 469   if(!bCell)
 470     Error("DigitizeEnergy","Wrong cell id number") ;
 471   geom->GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta);
 472
 473   if(emcalLoader->CalibData()) {
 474     fADCpedestalEC = emcalLoader->CalibData()
 475       ->GetADCpedestal(iSupMod,ieta,iphi);
 476     fADCchannelEC = emcalLoader->CalibData()
 477       ->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   else if (fManager)
 520     fLastEvent = fFirstEvent ; // what is this ??
 521
 522   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 523   Int_t ievent;
 524
 525   rl->LoadSDigits("EMCAL");
 526   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 527
 528     rl->GetEvent(ievent);
 529
 530     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 531
 532     WriteDigits() ;
 533
 534     if(strstr(option,"deb"))
 535       PrintDigits(option);
 536     if(strstr(option,"table")) gObjectTable->Print();
 537
 538     //increment the total number of Digits per run
 539     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
 540   }
 541
 542   emcalLoader->CleanDigitizer() ;
 543
 544   if(strstr(option,"tim")){
 545     gBenchmark->Stop("EMCALDigitizer");
 546     AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event",
 547          gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
 548   }
 549 }
 550
 551 //____________________________________________________________________________
 552 Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks)
 553 {
 554   //  Returns the shortest time among all time ticks
 555
 556   ticks->Sort() ; //Sort in accordance with times of ticks
 557   TIter it(ticks) ;
 558   AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
 559   Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 560
 561   AliEMCALTick * t ;
 562   while((t=(AliEMCALTick*) it.Next())){
 563     if(t->GetTime() < time)  //This tick starts before crossing
 564       *ctick+=*t ;
 565     else
 566       return time ;
 567
 568     time = ctick->CrossingTime(fTimeThreshold) ;
 569   }
 570   return time ;
 571 }
 572
 573 //____________________________________________________________________________
 574 Bool_t AliEMCALDigitizer::Init()
 575 {
 576   // Makes all memory allocations
 577   fInit = kTRUE ;
 578   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 579
 580   if ( emcalLoader == 0 ) {
 581     Fatal("Init", "Could not obtain the AliEMCALLoader");
 582     return kFALSE;
 583   }
 584
 585   fFirstEvent = 0 ;
 586   fLastEvent = fFirstEvent ;
 587
 588   if(fManager)
 589     fInput = fManager->GetNinputs() ;
 590   else
 591     fInput           = 1 ;
 592
 593   fInputFileNames  = new TString[fInput] ;
 594   fEventNames      = new TString[fInput] ;
 595   fInputFileNames[0] = GetTitle() ;
 596   fEventNames[0]     = fEventFolderName.Data() ;
 597   Int_t index ;
 598   for (index = 1 ; index < fInput ; index++) {
 599     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 600     TString tempo = fManager->GetInputFolderName(index) ;
 601     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager
 602   }
 603
 604   //to prevent cleaning of this object while GetEvent is called
 605   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 606
 607   //PH  Print();
 608
 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 }