PHOS/AliPHOSDigitizer.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 /* History of cvs commits:
  19  *
  20  * $Log$
  21  * Revision 1.101  2007/10/14 21:08:10  schutz
  22  * Introduced the checking of QA results from previous step before entering the event loop
  23  *
  24  * Revision 1.100  2007/10/10 09:05:10  schutz
  25  * Changing name QualAss to QA
  26  *
  27  * Revision 1.99  2007/09/30 17:08:20  schutz
  28  * Introducing the notion of QA data acquisition cycle (needed by online)
  29  *
  30  * Revision 1.98  2007/09/26 14:22:17  cvetan
  31  * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
  32  *
  33  * Revision 1.97  2007/08/07 14:12:03  kharlov
  34  * Quality assurance added (Yves Schutz)
  35  *
  36  * Revision 1.96  2007/04/28 10:43:36  policheh
  37  * Dead channels simulation: digit energy sets to 0.
  38  *
  39  * Revision 1.95  2007/04/10 07:20:52  kharlov
  40  * Decalibration should use the same CDB as calibration in AliPHOSClusterizerv1
  41  *
  42  * Revision 1.94  2007/02/01 10:34:47  hristov
  43  * Removing warnings on Solaris x86
  44  *
  45  * Revision 1.93  2006/10/17 13:17:01  kharlov
  46  * Replace AliInfo by AliDebug
  47  *
  48  * Revision 1.92  2006/08/28 10:01:56  kharlov
  49  * Effective C++ warnings fixed (Timur Pocheptsov)
  50  *
  51  * Revision 1.91  2006/04/29 20:25:30  hristov
  52  * Decalibration is implemented (Yu.Kharlov)
  53  *
  54  * Revision 1.90  2006/04/22 10:30:17  hristov
  55  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
  56  *
  57  * Revision 1.89  2006/04/11 15:22:59  hristov
  58  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
  59  *
  60  * Revision 1.88  2006/03/13 14:05:43  kharlov
  61  * Calibration objects for EMC and CPV
  62  *
  63  * Revision 1.87  2005/08/24 15:33:49  kharlov
  64  * Calibration data for raw digits
  65  *
  66  * Revision 1.86  2005/07/12 20:07:35  hristov
  67  * Changes needed to run simulation and reconstrruction in the same AliRoot session
  68  *
  69  * Revision 1.85  2005/05/28 14:19:04  schutz
  70  * Compilation warnings fixed by T.P.
  71  *
  72  */
  73
  74 //_________________________________________________________________________
  75 //*-- Author :  Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  76 //////////////////////////////////////////////////////////////////////////////
  77 // This TTask performs digitization of Summable digits (in the PHOS case it is just
  78 // the sum of contributions from all primary particles into a given cell).
  79 // In addition it performs mixing of summable digits from different events.
  80 // The name of the TTask is also the title of the branch that will contain
  81 // the created SDigits
  82 // The title of the TTAsk is the name of the file that contains the hits from
  83 // which the SDigits are created
  84 //
  85 // For each event two branches are created in TreeD:
  86 //   "PHOS" - list of digits
  87 //   "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
  88 //
  89 // Note, that one can set a title for new digits branch, and repeat digitization with
  90 // another set of parameters.
  91 //
  92 // Use case:
  93 // root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
  94 // root[1] d->ExecuteTask()
  95 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  96 //                       //Digitizes SDigitis in all events found in file galice.root
  97 //
  98 // root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;
  99 //                       // Will read sdigits from galice1.root
 100 // root[3] d1->MixWith("galice2.root")
 101 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
 102 //                       // Reads another set of sdigits from galice2.root
 103 // root[3] d1->MixWith("galice3.root")
 104 //                       // Reads another set of sdigits from galice3.root
 105 // root[4] d->ExecuteTask("deb timing")
 106 //                       // Reads SDigits from files galice1.root, galice2.root ....
 107 //                       // mixes them and stores produced Digits in file galice1.root
 108 //                       // deb - prints number of produced digits
 109 //                       // deb all - prints list of produced digits
 110 //                       // timing  - prints time used for digitization
 111 //
 112
 113 // --- ROOT system ---
 114 #include "TTree.h"
 115 #include "TSystem.h"
 116 #include "TBenchmark.h"
 117 #include "TRandom.h"
 118
 119 // --- Standard library ---
 120
 121 // --- AliRoot header files ---
 122 #include "AliLog.h"
 123 #include "AliRunDigitizer.h"
 124 #include "AliPHOSDigit.h"
 125 #include "AliPHOSGetter.h"
 126 #include "AliPHOSDigitizer.h"
 127 #include "AliPHOSSDigitizer.h"
 128 #include "AliPHOSGeometry.h"
 129 #include "AliPHOSTick.h"
 130 #include "AliPHOSQADataMaker.h"
 131
 132 ClassImp(AliPHOSDigitizer)
 133
 134
 135 //____________________________________________________________________________
 136 AliPHOSDigitizer::AliPHOSDigitizer() :
 137   AliDigitizer("",""),
 138   fDefaultInit(kTRUE),
 139   fDigitsInRun(0),
 140   fInit(kFALSE),
 141   fInput(0),
 142   fInputFileNames(0x0),
 143   fEventNames(0x0),
 144   fEmcCrystals(0),
 145   fPinNoise(0.f),
 146   fEMCDigitThreshold(0.f),
 147   fCPVNoise(0.f),
 148   fCPVDigitThreshold(0.f),
 149   fTimeResolution(0.f),
 150   fTimeThreshold(0.f),
 151   fTimeSignalLength(0.f),
 152   fADCchanelEmc(0.f),
 153   fADCpedestalEmc(0.f),
 154   fNADCemc(0),
 155   fADCchanelCpv(0.f),
 156   fADCpedestalCpv(0.f),
 157   fNADCcpv(0),
 158   fEventFolderName(""),
 159   fFirstEvent(0),
 160   fLastEvent(0),
 161 //  fQADM (0x0),
 162   fEventCounter(0)
 163 {
 164   // ctor
 165   InitParameters() ;
 166   fManager = 0 ;                     // We work in the standalong mode
 167 }
 168
 169 //____________________________________________________________________________
 170 AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName,
 171                                    TString eventFolderName):
 172   AliDigitizer("PHOS"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
 173   fDefaultInit(kFALSE),
 174   fDigitsInRun(0),
 175   fInit(kFALSE),
 176   fInput(0),
 177   fInputFileNames(0x0),
 178   fEventNames(0x0),
 179   fEmcCrystals(0),
 180   fPinNoise(0.f),
 181   fEMCDigitThreshold(0.f),
 182   fCPVNoise(0.f),
 183   fCPVDigitThreshold(0.f),
 184   fTimeResolution(0.f),
 185   fTimeThreshold(0.f),
 186   fTimeSignalLength(0.f),
 187   fADCchanelEmc(0.f),
 188   fADCpedestalEmc(0.f),
 189   fNADCemc(0),
 190   fADCchanelCpv(0.f),
 191   fADCpedestalCpv(0.f),
 192   fNADCcpv(0),
 193   fEventFolderName(eventFolderName),
 194   fFirstEvent(0),
 195   fLastEvent(0),
 196 //  fQADM (0x0),
 197   fEventCounter(0)
 198 {
 199   // ctor
 200   InitParameters() ;
 201   Init() ;
 202   fDefaultInit = kFALSE ;
 203   fManager = 0 ;                     // We work in the standalong mode
 204  // //Initialize the quality assurance data maker only once
 205 //  fQADM = new AliPHOSQADataMaker() ;
 206 //  //FIXME: get the run number
 207 //  Int_t run = 0 ;
 208 //  //EMXIF
 209 //   GetQADataMaker()->Init(AliQA::kDIGITS, run, fgkCycles) ;
 210 //   GetQADataMaker()->StartOfCycle(AliQA::kDIGITS) ;
 211 }
 212
 213 //____________________________________________________________________________
 214 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) :
 215   AliDigitizer(d),
 216   fDefaultInit(d.fDefaultInit),
 217   fDigitsInRun(d.fDigitsInRun),
 218   fInit(d.fInit),
 219   fInput(d.fInput),
 220   fInputFileNames(0x0),//?
 221   fEventNames(0x0),//?
 222   fEmcCrystals(d.fEmcCrystals),
 223   fPinNoise(d.fPinNoise),
 224   fEMCDigitThreshold(d.fEMCDigitThreshold),
 225   fCPVNoise(d.fCPVNoise),
 226   fCPVDigitThreshold(d.fCPVDigitThreshold),
 227   fTimeResolution(d.fTimeResolution),
 228   fTimeThreshold(d.fTimeThreshold),
 229   fTimeSignalLength(d.fTimeSignalLength),
 230   fADCchanelEmc(d.fADCchanelEmc),
 231   fADCpedestalEmc(d.fADCpedestalEmc),
 232   fNADCemc(d.fNADCemc),
 233   fADCchanelCpv(d.fADCchanelCpv),
 234   fADCpedestalCpv(d.fADCpedestalCpv),
 235   fNADCcpv(d.fNADCcpv),
 236   fEventFolderName(d.fEventFolderName),
 237   fFirstEvent(d.fFirstEvent),
 238   fLastEvent(d.fLastEvent),
 239 //  fQADM (d.fQADM),
 240   fEventCounter(0)
 241
 242 {
 243   // copyy ctor
 244   SetName(d.GetName()) ;
 245   SetTitle(d.GetTitle()) ;
 246 //Initialize the quality assurance data maker only once
 247  //FIXME: get the run number
 248 //  Int_t run = 0 ;
 249 //  //EMXIF
 250 //  GetQADataMaker()->Init(AliQA::kDIGITS, run, fgkCycles) ;
 251 //  GetQADataMaker()->StartOfCycle(AliQA::kDIGITS) ;
 252 }
 253
 254 //____________________________________________________________________________
 255 AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd) :
 256   AliDigitizer(rd,"PHOS"+AliConfig::Instance()->GetDigitizerTaskName()),
 257   fDefaultInit(kFALSE),
 258   fDigitsInRun(0),
 259   fInit(kFALSE),
 260   fInput(0),
 261   fInputFileNames(0x0),
 262   fEventNames(0x0),
 263   fEmcCrystals(0),
 264   fPinNoise(0.f),
 265   fEMCDigitThreshold(0.f),
 266   fCPVNoise(0.f),
 267   fCPVDigitThreshold(0.f),
 268   fTimeResolution(0.f),
 269   fTimeThreshold(0.f),
 270   fTimeSignalLength(0.f),
 271   fADCchanelEmc(0.f),
 272   fADCpedestalEmc(0.f),
 273   fNADCemc(0),
 274   fADCchanelCpv(0.f),
 275   fADCpedestalCpv(0.f),
 276   fNADCcpv(0),
 277   fEventFolderName(fManager->GetInputFolderName(0)),
 278   fFirstEvent(0),
 279   fLastEvent(0),
 280 //  fQADM (0x0),
 281   fEventCounter(0)
 282
 283 {
 284   // ctor Init() is called by RunDigitizer
 285   fManager = rd ;
 286   SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
 287   InitParameters() ;
 288   fDefaultInit = kFALSE ;
 289 //Initialize the quality assurance data maker only once
 290 //  fQADM = new AliPHOSQADataMaker() ;
 291 // //FIXME: get the run number
 292 //  Int_t run = 0 ;
 293 //  //EMXIF
 294 //  GetQADataMaker()->Init(AliQA::kDIGITS, run) ;
 295 //  GetQADataMaker()->StartOfCycle(AliQA::kDIGITS) ;
 296 }
 297
 298 //____________________________________________________________________________
 299   AliPHOSDigitizer::~AliPHOSDigitizer()
 300 {
 301   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ;
 302
 303   // Clean Digitizer from the white board
 304   gime->PhosLoader()->CleanDigitizer() ;
 305   // dtor
 306   delete [] fInputFileNames ;
 307   delete [] fEventNames ;
 308
 309 //  delete fQADM ;
 310
 311 }
 312
 313 //____________________________________________________________________________
 314 void AliPHOSDigitizer::Digitize(Int_t event)
 315 {
 316
 317   // Makes the digitization of the collected summable digits.
 318   //  It first creates the array of all PHOS modules
 319   //  filled with noise (different for EMC, and CPV) and
 320   //  then adds contributions from SDigits.
 321   // This design avoids scanning over the list of digits to add
 322   // contribution to new SDigits only.
 323
 324   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ;
 325   Int_t ReadEvent = event ;
 326   if (fManager)
 327     ReadEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
 328   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 329                   ReadEvent, GetTitle(), fEventFolderName.Data())) ;
 330   gime->Event(ReadEvent, "S") ;
 331   TClonesArray * digits = gime->Digits() ;
 332   digits->Clear() ;
 333
 334   const AliPHOSGeometry *geom = gime->PHOSGeometry() ;
 335   //Making digits with noise, first EMC
 336   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 337
 338   Int_t nCPV ;
 339   Int_t absID ;
 340
 341   nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * geom->GetNModules() ;
 342
 343   digits->Expand(nCPV) ;
 344
 345   // get first the sdigitizer from the tasks list
 346   if ( !gime->SDigitizer() )
 347     gime->LoadSDigitizer();
 348   AliPHOSSDigitizer * sDigitizer = gime->SDigitizer();
 349
 350   if ( !sDigitizer )
 351     AliFatal(Form("SDigitizer with name %s %s not found",
 352                   GetTitle(), fEventFolderName.Data() )) ;
 353
 354   //take all the inputs to add together and load the SDigits
 355   TObjArray * sdigArray = new TObjArray(fInput) ;
 356   sdigArray->AddAt(gime->SDigits(), 0) ;
 357   Int_t i ;
 358   for(i = 1 ; i < fInput ; i++){
 359     TString tempo(fEventNames[i]) ;
 360     tempo += i ;
 361     AliPHOSGetter * gime1 = AliPHOSGetter::Instance(fInputFileNames[i], tempo) ;
 362     if (fManager)
 363       ReadEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
 364     AliInfo(Form("Adding event %d from input stream %d %s %s",
 365                  ReadEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
 366     gime1->Event(ReadEvent,"S");
 367     sdigArray->AddAt(gime1->SDigits(), i) ;
 368   }
 369
 370   //Find the first crystal with signal
 371   Int_t nextSig = 200000 ;
 372   TClonesArray * sdigits ;
 373   for(i = 0 ; i < fInput ; i++){
 374     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 375     if ( !sdigits->GetEntriesFast() )
 376       continue ;
 377     Int_t curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
 378     if(curNext < nextSig)
 379       nextSig = curNext ;
 380   }
 381
 382   TArrayI index(fInput) ;
 383   index.Reset() ;  //Set all indexes to zero
 384
 385   AliPHOSDigit * digit ;
 386   AliPHOSDigit * curSDigit ;
 387
 388   TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
 389
 390   //Put Noise contribution
 391   for(absID = 1 ; absID <= nEMC ; absID++){
 392     Float_t noise = gRandom->Gaus(0., fPinNoise) ;
 393     // YVK: do not digitize amplitudes for EMC
 394 //     new((*digits)[absID-1]) AliPHOSDigit( -1, absID, sDigitizer->Digitize(noise), TimeOfNoise() ) ;
 395     new((*digits)[absID-1]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
 396     //look if we have to add signal?
 397     digit = dynamic_cast<AliPHOSDigit *>(digits->At(absID-1)) ;
 398
 399     if(absID==nextSig){
 400       //Add SDigits from all inputs
 401       ticks->Clear() ;
 402       Int_t contrib = 0 ;
 403       Float_t a = digit->GetEnergy() ;
 404       Float_t b = TMath::Abs( a / fTimeSignalLength) ;
 405       //Mark the beginning of the signal
 406       new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);
 407       //Mark the end of the signal
 408       new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b);
 409
 410       //loop over inputs
 411       for(i = 0 ; i < fInput ; i++){
 412         if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 413           curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 414         else
 415           curSDigit = 0 ;
 416         //May be several digits will contribute from the same input
 417         while(curSDigit && curSDigit->GetId() == absID){
 418           //Shift primary to separate primaries belonging different inputs
 419           Int_t primaryoffset ;
 420           if(fManager)
 421             primaryoffset = fManager->GetMask(i) ;
 422           else
 423             primaryoffset = 10000000*i ;
 424           curSDigit->ShiftPrimary(primaryoffset) ;
 425
 426           a = curSDigit->GetEnergy() ;
 427           b = a /fTimeSignalLength ;
 428           new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);
 429           new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
 430
 431           *digit += *curSDigit ;  //add energies
 432
 433           index[i]++ ;
 434           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 435             curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 436           else
 437             curSDigit = 0 ;
 438         }
 439       }
 440
 441       //calculate and set time
 442       Float_t time = FrontEdgeTime(ticks) ;
 443       digit->SetTime(time) ;
 444
 445       //Find next signal module
 446       nextSig = 200000 ;
 447       for(i = 0 ; i < fInput ; i++){
 448         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 449         Int_t curNext = nextSig ;
 450         if(sdigits->GetEntriesFast() > index[i] ){
 451           curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
 452         }
 453         if(curNext < nextSig) nextSig = curNext ;
 454       }
 455     }
 456   }
 457
 458   ticks->Delete() ;
 459   delete ticks ;
 460
 461   //Now CPV digits (different noise and no timing)
 462   for(absID = nEMC+1; absID <= nCPV; absID++){
 463     Float_t noise = gRandom->Gaus(0., fCPVNoise) ;
 464     new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
 465     //look if we have to add signal?
 466     if(absID==nextSig){
 467       digit = dynamic_cast<AliPHOSDigit *>(digits->At(absID-1)) ;
 468       //Add SDigits from all inputs
 469       for(i = 0 ; i < fInput ; i++){
 470         if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 471           curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 472         else
 473           curSDigit = 0 ;
 474
 475         //May be several digits will contribute from the same input
 476         while(curSDigit && curSDigit->GetId() == absID){
 477           //Shift primary to separate primaries belonging different inputs
 478           Int_t primaryoffset ;
 479           if(fManager)
 480             primaryoffset = fManager->GetMask(i) ;
 481           else
 482             primaryoffset = 10000000*i ;
 483           curSDigit->ShiftPrimary(primaryoffset) ;
 484
 485           //add energies
 486           *digit += *curSDigit ;
 487           index[i]++ ;
 488           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 489             curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;
 490           else
 491             curSDigit = 0 ;
 492         }
 493       }
 494
 495       //Find next signal module
 496       nextSig = 200000 ;
 497       for(i = 0 ; i < fInput ; i++){
 498         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 499         Int_t curNext = nextSig ;
 500         if(sdigits->GetEntriesFast() > index[i] )
 501           curNext = dynamic_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
 502         if(curNext < nextSig) nextSig = curNext ;
 503       }
 504
 505     }
 506   }
 507
 508   delete sdigArray ; //We should not delete its contents
 509
 510   //remove digits below thresholds
 511   for(i = 0 ; i < nEMC ; i++){
 512     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
 513     DecalibrateEMC(digit);
 514     if(digit->GetEnergy() < fEMCDigitThreshold)
 515       digits->RemoveAt(i) ;
 516     else
 517       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
 518   }
 519
 520
 521   for(i = nEMC; i < nCPV ; i++)
 522 //     if( sDigitizer->Calibrate( dynamic_cast<AliPHOSDigit*>(digits->At(i))->GetAmp() ) < fCPVDigitThreshold )
 523     if( dynamic_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < fCPVDigitThreshold )
 524       digits->RemoveAt(i) ;
 525
 526   digits->Compress() ;
 527
 528   Int_t ndigits = digits->GetEntriesFast() ;
 529   digits->Expand(ndigits) ;
 530
 531   //Set indexes in list of digits and make true digitization of the energy
 532   for (i = 0 ; i < ndigits ; i++) {
 533     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
 534     digit->SetIndexInList(i) ;
 535     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
 536       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
 537     }
 538   }
 539
 540   Int_t relId[4];
 541
 542   //set amplitudes in bad channels to zero
 543   for(i = 0 ; i <digits->GetEntries(); i++){
 544     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
 545     gime->PHOSGeometry()->AbsToRelNumbering(digit->GetId(),relId);
 546     if(relId[1] == 0) // Emc
 547       if(gime->CalibData()->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.);
 548   }
 549
 550 }
 551
 552 //____________________________________________________________________________
 553 void AliPHOSDigitizer::DecalibrateEMC(AliPHOSDigit *digit)
 554 {
 555   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
 556
 557   AliPHOSGetter* gime = AliPHOSGetter::Instance();
 558
 559   if(!gime->CalibData()) {
 560     AliPHOSCalibData* cdb = new AliPHOSCalibData(-1);
 561     gime->SetCalibData(cdb);
 562   }
 563
 564   //Determine rel.position of the cell absolute ID
 565   Int_t relId[4];
 566   gime->PHOSGeometry()->AbsToRelNumbering(digit->GetId(),relId);
 567   Int_t module=relId[0];
 568   Int_t row   =relId[2];
 569   Int_t column=relId[3];
 570   Float_t decalibration = gime->CalibData()->GetADCchannelEmc(module,column,row);
 571   Float_t energy = digit->GetEnergy() / decalibration;
 572   digit->SetEnergy(energy);
 573 }
 574 //____________________________________________________________________________
 575 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
 576 {
 577   // Returns digitized value of the CPV charge in a pad absId
 578
 579   AliPHOSGetter* gime = AliPHOSGetter::Instance();
 580
 581   if(!gime->CalibData()) {
 582     AliPHOSCalibData* cdb = new AliPHOSCalibData(-1); // use AliCDBManager's run number
 583     gime->SetCalibData(cdb);
 584   }
 585
 586   //Determine rel.position of the cell absId
 587   Int_t relId[4];
 588   gime->PHOSGeometry()->AbsToRelNumbering(absId,relId);
 589   Int_t module=relId[0];
 590   Int_t row   =relId[2];
 591   Int_t column=relId[3];
 592
 593   Int_t channel = 0;
 594
 595   if(absId > fEmcCrystals){ //digitize CPV only
 596
 597     //reading calibration data for cell absId.
 598     //If no calibration DB found, accept default values.
 599
 600     if(gime->CalibData()) {
 601       fADCpedestalCpv = gime->CalibData()->GetADCpedestalCpv(module,column,row);
 602       fADCchanelCpv   = gime->CalibData()->GetADCchannelCpv( module,column,row);
 603     }
 604
 605     channel = (Int_t) TMath::Ceil((charge - fADCpedestalCpv)/fADCchanelCpv) ;
 606     if(channel > fNADCcpv ) channel =  fNADCcpv ;
 607   }
 608   return channel ;
 609 }
 610
 611 //____________________________________________________________________________
 612 void AliPHOSDigitizer::Exec(Option_t *option)
 613 {
 614   // Steering method to process digitization for events
 615   // in the range from fFirstEvent to fLastEvent.
 616   // This range is optionally set by SetEventRange().
 617   // if fLastEvent=-1, then process events until the end.
 618   // by default fLastEvent = fFirstEvent (process only one event)
 619
 620   if (!fInit) { // to prevent overwrite existing file
 621     AliError(Form("Give a version name different from %s",
 622                   fEventFolderName.Data() )) ;
 623     return ;
 624   }
 625
 626   if (strstr(option,"print")) {
 627     Print();
 628     return ;
 629   }
 630
 631   // check the QA result for Hits and SDigits
 632   AliQA * qa = AliQA::Instance(AliQA::kPHOS) ;
 633   if ( qa->IsSet(AliQA::kPHOS, AliQA::kSIM, AliQA::kFATAL)) {
 634         AliFatal("QA status in Hits and/or SDIGITS was Fatal") ;
 635   } else if ( qa->IsSet(AliQA::kPHOS, AliQA::kSIM, AliQA::kERROR)) {
 636         AliError("QA status in Hits and/or SDIGITS was Error") ;
 637   } else if ( qa->IsSet(AliQA::kPHOS, AliQA::kSIM, AliQA::kWARNING) ) {
 638         AliWarning("QA status in Hits and/or SDIGITS was Warning") ;
 639   } else if ( qa->IsSet(AliQA::kPHOS, AliQA::kSIM, AliQA::kINFO) ) {
 640         AliInfo("QA status in Hits and/or SDIGITS was Info") ;
 641   }
 642
 643  if(strstr(option,"tim"))
 644      gBenchmark->Start("PHOSDigitizer");
 645
 646   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ;
 647
 648   // Post Digitizer to the white board
 649   gime->PostDigitizer(this) ;
 650
 651   if (fLastEvent == -1)
 652     fLastEvent = gime->MaxEvent() - 1 ;
 653   else if (fManager)
 654     fLastEvent = fFirstEvent ;
 655
 656   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 657
 658   Int_t ievent ;
 659
 660   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 661     fEventCounter++ ;
 662     gime->Event(ievent,"S") ;
 663
 664     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 665
 666 //    //makes the quality assurance data
 667 //    if (GetQADataMaker()->IsCycleDone() ) {
 668 //      GetQADataMaker()->EndOfCycle(AliQA::kDIGITS) ;
 669 //        GetQADataMaker()->StartOfCycle(AliQA::kDIGITS) ;
 670 //    }
 671 //    GetQADataMaker()->Exec(AliQA::kDIGITS, gime->Digits()) ;
 672 //      GetQADataMaker()->Increment() ;
 673
 674     WriteDigits() ;
 675
 676     if(strstr(option,"deb"))
 677       PrintDigits(option);
 678
 679     //increment the total number of Digits per run
 680     fDigitsInRun += gime->Digits()->GetEntriesFast() ;
 681  }
 682
 683 //   //Write the quality assurance data only after the last event
 684 // if ( fEventCounter == gime->MaxEvent() ) {
 685 //    GetQADataMaker()->EndOfCycle(AliQA::kDIGITS) ;
 686 //      GetQADataMaker()->Finish() ;
 687 //  }
 688
 689   gime->PhosLoader()->CleanDigitizer();
 690
 691   if(strstr(option,"tim")){
 692     gBenchmark->Stop("PHOSDigitizer");
 693     TString message ;
 694     message = "  took %f seconds for Digitizing %f seconds per event\n" ;
 695     AliInfo(Form( message.Data(),
 696          gBenchmark->GetCpuTime("PHOSDigitizer"),
 697          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
 698   }
 699 }
 700
 701 //____________________________________________________________________________
 702 Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
 703 {
 704   // Returns the shortest time among all time ticks
 705
 706   ticks->Sort() ; //Sort in accordance with times of ticks
 707   TIter it(ticks) ;
 708   AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
 709   Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 710
 711   AliPHOSTick * t ;
 712   while((t=(AliPHOSTick*) it.Next())){
 713     if(t->GetTime() < time)  //This tick starts before crossing
 714       *ctick+=*t ;
 715     else
 716       return time ;
 717
 718     time = ctick->CrossingTime(fTimeThreshold) ;
 719   }
 720   return time ;
 721 }
 722
 723 //____________________________________________________________________________
 724 Bool_t AliPHOSDigitizer::Init()
 725 {
 726   // Makes all memory allocations
 727   fInit = kTRUE ;
 728   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 729   if ( gime == 0 ) {
 730     AliFatal(Form("Could not obtain the Getter object for file %s and event %s !",
 731                   GetTitle(), fEventFolderName.Data()));
 732     return kFALSE;
 733   }
 734
 735   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 736
 737   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
 738
 739   TString opt("Digits") ;
 740   if(gime->VersionExists(opt) ) {
 741     AliError(Form("Give a version name different from %s",
 742                   fEventFolderName.Data() )) ;
 743     fInit = kFALSE ;
 744   }
 745
 746   fFirstEvent = 0 ;
 747   fLastEvent = fFirstEvent ;
 748   if (fManager)
 749     fInput = fManager->GetNinputs() ;
 750   else
 751     fInput           = 1 ;
 752
 753   fInputFileNames  = new TString[fInput] ;
 754   fEventNames      = new TString[fInput] ;
 755   fInputFileNames[0] = GetTitle() ;
 756   fEventNames[0]     = fEventFolderName.Data() ;
 757   Int_t index ;
 758   for (index = 1 ; index < fInput ; index++) {
 759     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 760     TString tempo = fManager->GetInputFolderName(index) ;
 761     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
 762   }
 763
 764   //to prevent cleaning of this object while GetEvent is called
 765   gime->PhosLoader()->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 766
 767   return fInit ;
 768 }
 769
 770 //____________________________________________________________________________
 771 void AliPHOSDigitizer::InitParameters()
 772 {
 773   // Set initial parameters Digitizer
 774
 775   fPinNoise           = 0.004 ;  // [GeV]
 776   fEMCDigitThreshold  = 0.012 ;  // [GeV]
 777   fCPVNoise           = 0.01;    // [aux units]
 778   fCPVDigitThreshold  = 0.09 ;   // [aux units]
 779   fTimeResolution     = 0.5e-9 ; // [sec]
 780   fTimeSignalLength   = 1.0e-9 ; // [sec]
 781   fDigitsInRun  = 0 ;
 782   fADCchanelEmc = 1.0;        // Coefficient between real and measured energies in EMC
 783   fADCpedestalEmc = 0. ;      //
 784   fNADCemc = (Int_t) TMath::Power(2,16) ;  // number of channels in EMC ADC
 785
 786   fADCchanelCpv = 0.0012 ;          // width of one ADC channel in CPV 'popugais'
 787   fADCpedestalCpv = 0.012 ;         //
 788   fNADCcpv = (Int_t) TMath::Power(2,12);      // number of channels in CPV ADC
 789
 790 //   fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
 791   fTimeThreshold = 0.001 ; // [GeV]
 792   SetEventRange(0,-1) ;
 793
 794 }
 795
 796 //__________________________________________________________________
 797 void AliPHOSDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
 798 {
 799   // Allows to produce digits by superimposing background and signal event.
 800   // It is assumed, that headers file with SIGNAL events is opened in
 801   // the constructor.
 802   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed
 803   // Thus we avoid writing (changing) huge and expensive
 804   // backgound files: all output will be writen into SIGNAL, i.e.
 805   // opened in constructor file.
 806   //
 807   // One can open as many files to mix with as one needs.
 808   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
 809   // can be mixed.
 810
 811   if( strcmp(fEventFolderName, "") == 0 )
 812     Init() ;
 813
 814   if(fManager){
 815     Warning("MixWith", "Cannot use this method with AliRunDigitizer\n" ) ;
 816     return ;
 817   }
 818   // looking for file which contains AliRun
 819   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
 820     AliError(Form("File %s does not exist!", alirunFileName.Data())) ;
 821     return ;
 822   }
 823   // looking for the file which contains SDigits
 824   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 825   TString fileName( gime->GetSDigitsFileName() ) ;
 826     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 827       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
 828     if ( (gSystem->AccessPathName(fileName)) ) {
 829       AliError(Form("The file %s does not exist!", fileName.Data())) ;
 830       return ;
 831     }
 832     // need to increase the arrays
 833     TString tempo = fInputFileNames[fInput-1] ;
 834     delete [] fInputFileNames ;
 835     fInputFileNames = new TString[fInput+1] ;
 836     fInputFileNames[fInput-1] = tempo ;
 837
 838     tempo = fEventNames[fInput-1] ;
 839     delete [] fEventNames ;
 840     fEventNames = new TString[fInput+1] ;
 841     fEventNames[fInput-1] = tempo ;
 842
 843     fInputFileNames[fInput] = alirunFileName ;
 844     fEventNames[fInput]     = eventFolderName ;
 845     fInput++ ;
 846 }
 847
 848 //__________________________________________________________________
 849 void AliPHOSDigitizer::Print(const Option_t *)const
 850 {
 851   // Print Digitizer's parameters
 852   AliInfo(Form("\n------------------- %s -------------", GetName() )) ;
 853   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 854     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 855
 856     Int_t nStreams ;
 857     if (fManager)
 858       nStreams =  GetNInputStreams() ;
 859     else
 860       nStreams = fInput ;
 861
 862     Int_t index = 0 ;
 863     for (index = 0 ; index < nStreams ; index++) {
 864       TString tempo(fEventNames[index]) ;
 865       tempo += index ;
 866       AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileNames[index], tempo) ;
 867       TString fileName( gime->GetSDigitsFileName() ) ;
 868       if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
 869         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
 870       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
 871     }
 872     AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 873     printf("\nWriting digits to %s", gime->GetDigitsFileName().Data()) ;
 874
 875     printf("\nWith following parameters:\n") ;
 876     printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ;
 877     printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;
 878     printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ;
 879     printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;
 880     printf(" ---------------------------------------------------\n") ;
 881   }
 882   else
 883     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
 884
 885 }
 886
 887 //__________________________________________________________________
 888  void AliPHOSDigitizer::PrintDigits(Option_t * option)
 889 {
 890   // Print a table of digits
 891
 892   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 893   TClonesArray * digits = gime->Digits() ;
 894
 895   AliInfo(Form("%d", digits->GetEntriesFast())) ;
 896   printf("\nevent %d", gAlice->GetEvNumber()) ;
 897   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;
 898
 899
 900   if(strstr(option,"all")||strstr(option,"EMC")){
 901     //loop over digits
 902     AliPHOSDigit * digit;
 903     printf("\nEMC digits (with primaries):\n")  ;
 904     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 905     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
 906     Int_t index ;
 907     for (index = 0 ; (index < digits->GetEntriesFast()) &&
 908            (dynamic_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
 909       digit = (AliPHOSDigit * )  digits->At(index) ;
 910       if(digit->GetNprimary() == 0)
 911         continue;
 912 //       printf("%6d  %8d    %6.5e %4d      %2d :",
 913 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
 914       printf("%6d  %.4f    %6.5e %4d      %2d :",
 915               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 916       Int_t iprimary;
 917       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 918         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 919       }
 920       printf("\n") ;
 921     }
 922   }
 923
 924   if(strstr(option,"all")||strstr(option,"CPV")){
 925
 926     //loop over CPV digits
 927     AliPHOSDigit * digit;
 928     printf("\nCPV digits (with primaries):\n")  ;
 929     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 930     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
 931     Int_t index ;
 932     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
 933       digit = (AliPHOSDigit * )  digits->At(index) ;
 934       if(digit->GetId() > maxEmc){
 935         printf("%6d  %8d    %4d      %2d :",
 936                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
 937         Int_t iprimary;
 938         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 939           printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 940         }
 941         printf("\n") ;
 942       }
 943     }
 944   }
 945
 946 }
 947
 948 //__________________________________________________________________
 949 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
 950 {  // Calculates the time signal generated by noise
 951   //PH  Info("TimeOfNoise", "Change me") ;
 952   return gRandom->Rndm() * 1.28E-5;
 953 }
 954
 955 //__________________________________________________________________
 956 void AliPHOSDigitizer::Unload()
 957 {
 958
 959   Int_t i ;
 960   for(i = 1 ; i < fInput ; i++){
 961     TString tempo(fEventNames[i]) ;
 962     tempo += i ;
 963     AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileNames[i], tempo) ;
 964     gime->PhosLoader()->UnloadSDigits() ;
 965   }
 966
 967   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 968   gime->PhosLoader()->UnloadDigits() ;
 969 }
 970
 971 //____________________________________________________________________________
 972 void AliPHOSDigitizer::WriteDigits()
 973 {
 974
 975   // Makes TreeD in the output file.
 976   // Check if branch already exists:
 977   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
 978   //      already existing branches.
 979   //   else creates branch with Digits, named "PHOS", title "...",
 980   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
 981   //      and names of files, from which digits are made.
 982
 983   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ;
 984   const TClonesArray * digits = gime->Digits() ;
 985   TTree * treeD = gime->TreeD();
 986
 987   // -- create Digits branch
 988   Int_t bufferSize = 32000 ;
 989   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
 990   digitsBranch->SetTitle(fEventFolderName);
 991   digitsBranch->Fill() ;
 992
 993   gime->WriteDigits("OVERWRITE");
 994   gime->WriteDigitizer("OVERWRITE");
 995
 996   Unload() ;
 997
 998 }