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