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