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