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