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