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