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