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