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 <TGeoManager.h>
 133 #include "AliLog.h"
 134 #include "AliRunDigitizer.h"
 135 #include "AliPHOSDigit.h"
 136 #include "AliPHOSDigitizer.h"
 137 #include "AliPHOSGeometry.h"
 138 #include "AliPHOSTick.h"
 139 #include "AliPHOSSimParam.h"
 140 #include "AliPHOSCalibData.h"
 141 #include "AliRunLoader.h"
 142 #include "AliPHOSLoader.h"
 143 #include "AliPHOSPulseGenerator.h"
 144
 145 ClassImp(AliPHOSDigitizer)
 146
 147
 148 //____________________________________________________________________________
 149 AliPHOSDigitizer::AliPHOSDigitizer() :
 150   AliDigitizer("",""),
 151   fDefaultInit(kTRUE),
 152   fDigitsInRun(0),
 153   fInit(kFALSE),
 154   fInput(0),
 155   fInputFileNames(0x0),
 156   fEventNames(0x0),
 157   fEmcCrystals(0),
 158   fEventFolderName(""),
 159   fFirstEvent(0),
 160   fLastEvent(0),
 161   fcdb(0x0),
 162   fEventCounter(0),
 163   fPulse(0),
 164   fADCValuesLG(0),
 165   fADCValuesHG(0)
 166 {
 167   // ctor
 168   InitParameters() ;
 169   fManager = 0 ;                     // We work in the standalong mode
 170 }
 171
 172 //____________________________________________________________________________
 173 AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName,
 174                                    TString eventFolderName):
 175   AliDigitizer("PHOS"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
 176   fDefaultInit(kFALSE),
 177   fDigitsInRun(0),
 178   fInit(kFALSE),
 179   fInput(0),
 180   fInputFileNames(0x0),
 181   fEventNames(0x0),
 182   fEmcCrystals(0),
 183   fEventFolderName(eventFolderName),
 184   fFirstEvent(0),
 185   fLastEvent(0),
 186   fcdb(0x0),
 187   fEventCounter(0),
 188   fPulse(0),
 189   fADCValuesLG(0),
 190   fADCValuesHG(0)
 191 {
 192   // ctor
 193   InitParameters() ;
 194   Init() ;
 195   fDefaultInit = kFALSE ;
 196   fManager = 0 ;                     // We work in the standalone mode
 197   fcdb = new AliPHOSCalibData(-1);
 198 }
 199
 200 //____________________________________________________________________________
 201 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) :
 202   AliDigitizer(d),
 203   fDefaultInit(d.fDefaultInit),
 204   fDigitsInRun(d.fDigitsInRun),
 205   fInit(d.fInit),
 206   fInput(d.fInput),
 207   fInputFileNames(0x0),//?
 208   fEventNames(0x0),//?
 209   fEmcCrystals(d.fEmcCrystals),
 210   fEventFolderName(d.fEventFolderName),
 211   fFirstEvent(d.fFirstEvent),
 212   fLastEvent(d.fLastEvent),
 213   fcdb (0x0),
 214   fEventCounter(0),
 215   fPulse(0),
 216   fADCValuesLG(0),
 217   fADCValuesHG(0)
 218 {
 219   // copyy ctor
 220   SetName(d.GetName()) ;
 221   SetTitle(d.GetTitle()) ;
 222   fcdb = new AliPHOSCalibData(-1);
 223 }
 224
 225 //____________________________________________________________________________
 226 AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd) :
 227   AliDigitizer(rd,"PHOS"+AliConfig::Instance()->GetDigitizerTaskName()),
 228   fDefaultInit(kFALSE),
 229   fDigitsInRun(0),
 230   fInit(kFALSE),
 231   fInput(0),
 232   fInputFileNames(0x0),
 233   fEventNames(0x0),
 234   fEmcCrystals(0),
 235   fEventFolderName(fManager->GetInputFolderName(0)),
 236   fFirstEvent(0),
 237   fLastEvent(0),
 238   fcdb (0x0),
 239   fEventCounter(0),
 240   fPulse(0),
 241   fADCValuesLG(0),
 242   fADCValuesHG(0)
 243
 244 {
 245   // ctor Init() is called by RunDigitizer
 246   fManager = rd ;
 247   SetTitle(static_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
 248   InitParameters() ;
 249   fDefaultInit = kFALSE ;
 250   fcdb = new AliPHOSCalibData(-1);
 251 }
 252
 253 //____________________________________________________________________________
 254   AliPHOSDigitizer::~AliPHOSDigitizer()
 255 {
 256   // dtor
 257   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 258   if(rl){
 259     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 260
 261     if(phosLoader)
 262       phosLoader->CleanDigitizer() ;
 263   }
 264
 265   delete [] fInputFileNames ;
 266   delete [] fEventNames ;
 267
 268   delete fPulse;
 269   delete [] fADCValuesLG;
 270   delete [] fADCValuesHG;
 271
 272   if(fcdb){ delete fcdb ; fcdb=0;}
 273
 274 }
 275
 276 //____________________________________________________________________________
 277 void AliPHOSDigitizer::Digitize(Int_t event)
 278 {
 279
 280   // Makes the digitization of the collected summable digits.
 281   //  It first creates the array of all PHOS modules
 282   //  filled with noise (different for EMC, and CPV) and
 283   //  then adds contributions from SDigits.
 284   // This design avoids scanning over the list of digits to add
 285   // contribution to new SDigits only.
 286
 287   Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
 288
 289   //First stream
 290   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 291   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 292
 293   Int_t readEvent = event ;
 294   if (fManager)
 295     readEvent = static_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
 296   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 297                   readEvent, GetTitle(), fEventFolderName.Data())) ;
 298   rl->GetEvent(readEvent) ;
 299   phosLoader->CleanSDigits() ;
 300   phosLoader->LoadSDigits("READ") ;
 301
 302   //Prepare Output
 303   TClonesArray * digits = phosLoader->Digits() ;
 304   if( !digits ) {
 305     phosLoader->MakeDigitsArray() ;
 306     digits = phosLoader->Digits() ;
 307   }
 308
 309   digits->Clear() ;
 310
 311   //
 312   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 313   //Making digits with noise, first EMC
 314   //Check which PHOS modules are present
 315   Bool_t isPresent[5] ;
 316   TString volpath ;
 317   Int_t nmod=0 ;
 318   for(Int_t i=0; i<5; i++){
 319     volpath = "/ALIC_1/PHOS_";
 320     volpath += i+1;
 321     if (gGeoManager->CheckPath(volpath.Data())) {
 322       isPresent[i]=1 ;
 323       nmod++ ;
 324     }
 325     else{
 326       isPresent[i]=0 ;
 327     }
 328   }
 329
 330   Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
 331
 332   Int_t nCPV ;
 333   Int_t absID ;
 334
 335   //check if CPV exists
 336   Bool_t isCPVpresent=0 ;
 337   for(Int_t i=1; i<=5 && !isCPVpresent; i++){
 338     volpath = "/ALIC_1/PHOS_";
 339     volpath += i;
 340     volpath += "/PCPV_1";
 341     if (gGeoManager->CheckPath(volpath.Data()))
 342       isCPVpresent=1 ;
 343   }
 344
 345   if(isCPVpresent){
 346     nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
 347   }
 348   else{
 349      nCPV = nEMC ;
 350   }
 351
 352   digits->Expand(nCPV) ;
 353
 354   //take all the inputs to add together and load the SDigits
 355   TObjArray * sdigArray = new TObjArray(fInput) ;
 356   sdigArray->AddAt(phosLoader->SDigits(), 0) ;
 357
 358   for(Int_t i = 1 ; i < fInput ; i++){
 359     TString tempo(fEventNames[i]) ;
 360     tempo += i ;
 361     AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
 362     if(!rl2){
 363       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
 364       if (!rl2) {
 365         Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
 366         return ;
 367       }
 368       rl2->LoadHeader();
 369     }
 370     AliPHOSLoader * phosLoader2 = static_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
 371
 372     if(fManager){
 373       readEvent = static_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
 374     }
 375     TClonesArray * digs ;
 376     if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
 377       AliInfo(Form("Adding event %d from input stream %d %s %s",
 378                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
 379        rl2->GetEvent(readEvent) ;
 380        phosLoader2->CleanDigits() ;
 381        phosLoader2->LoadDigits("READ") ;
 382        digs = phosLoader2->Digits() ;
 383        toMakeNoise=0 ; //Do not add noise, it is already in stream
 384     }
 385     else{
 386       AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
 387                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
 388        rl2->GetEvent(readEvent) ;
 389        phosLoader2->CleanSDigits() ;
 390        phosLoader2->LoadSDigits("READ") ;
 391        digs = phosLoader2->SDigits() ;
 392     }
 393     sdigArray->AddAt(digs, i) ;
 394   }
 395
 396   //Find the first crystal with signal
 397   Int_t nextSig = 200000 ;
 398   TClonesArray * sdigits ;
 399   for(Int_t i = 0 ; i < fInput ; i++){
 400     sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 401     if ( !sdigits->GetEntriesFast() )
 402       continue ;
 403     Int_t curNext = static_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
 404     if(curNext < nextSig)
 405       nextSig = curNext ;
 406   }
 407
 408   TArrayI index(fInput) ;
 409   index.Reset() ;  //Set all indexes to zero
 410
 411   AliPHOSDigit * digit ;
 412   AliPHOSDigit * curSDigit ;
 413
 414 //  TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
 415
 416   //Put Noise contribution
 417   Double_t apdNoise = 0. ;
 418   if(toMakeNoise)
 419      apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ;
 420
 421   Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
 422   Int_t idigit= 0;
 423   for(Int_t imod=0; imod<5; imod++){
 424     if(!isPresent[imod])
 425       continue ;
 426     Int_t firstAbsId=imod*emcpermod+1 ;
 427     Int_t lastAbsId =(imod+1)*emcpermod ;
 428     for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
 429       Float_t noise = gRandom->Gaus(0.,apdNoise) ;
 430       new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
 431       //look if we have to add signal?
 432       digit = static_cast<AliPHOSDigit *>(digits->At(idigit)) ;
 433       idigit++ ;
 434
 435       if(absID==nextSig){
 436       //Add SDigits from all inputs
 437 //      ticks->Clear() ;
 438 //      Int_t contrib = 0 ;
 439
 440 //New Timing model is necessary
 441 //      Float_t a = digit->GetEnergy() ;
 442 //      Float_t b = TMath::Abs( a / fTimeSignalLength) ;
 443 //      //Mark the beginning of the signal
 444 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);
 445 //      //Mark the end of the signal
 446 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b);
 447
 448 // Calculate time as time of the largest digit
 449         Float_t time = digit->GetTime() ;
 450         Float_t eTime= digit->GetEnergy() ;
 451
 452         //loop over inputs
 453         for(Int_t i = 0 ; i < fInput ; i++){
 454           if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
 455             curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 456             if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
 457               curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
 458               curSDigit->SetTime(CalibrateT(curSDigit->GetTime(),curSDigit->GetId())) ;
 459             }
 460           }
 461           else
 462             curSDigit = 0 ;
 463           //May be several digits will contribute from the same input
 464           while(curSDigit && curSDigit->GetId() == absID){
 465             //Shift primary to separate primaries belonging different inputs
 466             Int_t primaryoffset ;
 467             if(fManager)
 468               primaryoffset = fManager->GetMask(i) ;
 469             else
 470               primaryoffset = 10000000*i ;
 471             curSDigit->ShiftPrimary(primaryoffset) ;
 472
 473 //New Timing model is necessary
 474 //        a = curSDigit->GetEnergy() ;
 475 //        b = a /fTimeSignalLength ;
 476 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);
 477 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
 478           if(curSDigit->GetEnergy()>eTime){
 479             eTime=curSDigit->GetEnergy() ;
 480             time=curSDigit->GetTime() ;
 481           }
 482             *digit += *curSDigit ;  //add energies
 483
 484             index[i]++ ;
 485             if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 486               curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 487             else
 488               curSDigit = 0 ;
 489           }
 490         }
 491
 492         //calculate and set time
 493 //New Timing model is necessary
 494 //      Float_t time = FrontEdgeTime(ticks) ;
 495         digit->SetTime(time) ;
 496
 497         //Find next signal module
 498         nextSig = 200000 ;
 499         for(Int_t i = 0 ; i < fInput ; i++){
 500           sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 501           Int_t curNext = nextSig ;
 502           if(sdigits->GetEntriesFast() > index[i] ){
 503             curNext = static_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
 504           }
 505           if(curNext < nextSig) nextSig = curNext ;
 506         }
 507       }
 508     }
 509   }
 510
 511
 512   //Apply non-linearity
 513   if(AliPHOSSimParam::GetInstance()->IsCellNonlinearityOn()){ //Apply non-lineairyt on cell level
 514     const Double_t aNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyA() ;
 515     const Double_t bNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyB() ;
 516     const Double_t cNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyC() ;
 517     for(Int_t i = 0 ; i < nEMC ; i++){
 518       digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 519       Double_t e= digit->GetEnergy() ;
 520       // version(1)      digit->SetEnergy(e*(1+a*TMath::Exp(-e/b))) ;
 521       digit->SetEnergy(e*cNL*(1.+aNL*TMath::Exp(-e*e/2./bNL/bNL))) ; //Better agreement with data...
 522     }
 523   }
 524
 525
 526   //distretize energy if necessary
 527   if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
 528     Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
 529     for(Int_t i = 0 ; i < nEMC ; i++){
 530       digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 531       digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
 532     }
 533   }
 534
 535   //Apply decalibration if necessary
 536   for(Int_t i = 0 ; i < nEMC ; i++){
 537     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 538     Decalibrate(digit) ;
 539   }
 540
 541 //  ticks->Delete() ;
 542 //  delete ticks ;
 543
 544   //Now CPV digits (different noise and no timing)
 545   Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
 546   Int_t nEMCtotal=emcpermod*5 ;
 547   Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
 548   if(isCPVpresent){  //CPV is present in current geometry
 549     for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
 550       if(!isPresent[imod])
 551         continue ;
 552       Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
 553       Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
 554       for(absID = firstAbsId; absID <= lastAbsId; absID++){
 555         Float_t noise = gRandom->Gaus(0., cpvNoise) ;
 556         new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
 557         idigit++ ;
 558         //look if we have to add signal?
 559         if(absID==nextSig){
 560           digit = static_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
 561           //Add SDigits from all inputs
 562           for(Int_t i = 0 ; i < fInput ; i++){
 563              if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 564                curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 565              else
 566                curSDigit = 0 ;
 567
 568              //May be several digits will contribute from the same input
 569              while(curSDigit && curSDigit->GetId() == absID){
 570                //Shift primary to separate primaries belonging different inputs
 571                Int_t primaryoffset ;
 572                if(fManager)
 573                  primaryoffset = fManager->GetMask(i) ;
 574                else
 575                  primaryoffset = 10000000*i ;
 576                curSDigit->ShiftPrimary(primaryoffset) ;
 577
 578                //add energies
 579                *digit += *curSDigit ;
 580                index[i]++ ;
 581                if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 582                  curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;
 583                else
 584                  curSDigit = 0 ;
 585              }
 586           }
 587
 588           //Find next signal module
 589           nextSig = 200000 ;
 590           for(Int_t i = 0 ; i < fInput ; i++){
 591             sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 592             Int_t curNext = nextSig ;
 593             if(sdigits->GetEntriesFast() > index[i] )
 594               curNext = static_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
 595             if(curNext < nextSig) nextSig = curNext ;
 596           }
 597
 598         }
 599       }
 600     }
 601   }
 602
 603   delete sdigArray ; //We should not delete its contents
 604
 605   Int_t relId[4];
 606
 607   //set amplitudes in bad channels to zero
 608
 609   for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
 610     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 611     geom->AbsToRelNumbering(digit->GetId(),relId);
 612     if(relId[1] == 0) // Emc
 613       if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.);
 614   }
 615
 616   //remove digits below thresholds
 617   Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
 618   for(Int_t i = 0 ; i < nEMC ; i++){
 619     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 620
 621     if(digit->GetEnergy() < emcThreshold){
 622       digits->RemoveAt(i) ;
 623       continue ;
 624     }
 625
 626     geom->AbsToRelNumbering(digit->GetId(),relId);
 627
 628 //    digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
 629
 630     Float_t tres = TimeResolution(digit->GetEnergy()) ;
 631     digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
 632
 633     fPulse->Reset();
 634     fPulse->SetAmplitude(digit->GetEnergy()/
 635                          fcdb->GetADCchannelEmc(relId[0],relId[3],relId[2]));
 636     fPulse->SetTZero(digit->GetTimeR());
 637     fPulse->MakeSamples();
 638     fPulse->GetSamples(fADCValuesHG, fADCValuesLG) ;
 639     Int_t nSamples = fPulse->GetRawFormatTimeBins();
 640     digit->SetALTROSamplesHG(nSamples,fADCValuesHG);
 641     digit->SetALTROSamplesLG(nSamples,fADCValuesLG);
 642   }
 643
 644   Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
 645   for(Int_t i = nEMC; i < nCPV ; i++){
 646     if( static_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
 647       digits->RemoveAt(i) ;
 648   }
 649
 650   digits->Compress() ;
 651   Int_t ndigits = digits->GetEntriesFast() ;
 652
 653   //Set indexes in list of digits and make true digitization of the energy
 654   for (Int_t i = 0 ; i < ndigits ; i++) {
 655     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 656     digit->SetIndexInList(i) ;
 657     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
 658       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
 659     }
 660   }
 661
 662 }
 663 //____________________________________________________________________________
 664 Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
 665   //Apply calibration
 666   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 667
 668   //Determine rel.position of the cell absolute ID
 669   Int_t relId[4];
 670   geom->AbsToRelNumbering(absId,relId);
 671   Int_t module=relId[0];
 672   Int_t row   =relId[2];
 673   Int_t column=relId[3];
 674   if(relId[1]==0){ //This Is EMC
 675     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
 676     return amp*calibration ;
 677   }
 678   return 0 ;
 679 }
 680 //____________________________________________________________________________
 681 void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
 682 {
 683   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
 684
 685   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 686
 687   //Determine rel.position of the cell absolute ID
 688   Int_t relId[4];
 689   geom->AbsToRelNumbering(digit->GetId(),relId);
 690   Int_t module=relId[0];
 691   Int_t row   =relId[2];
 692   Int_t column=relId[3];
 693   if(relId[1]==0){ //This Is EMC
 694     Float_t decalib     = fcdb->GetADCchannelEmcDecalib(module,column,row); // O(1)
 695     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row)*decalib;
 696     Float_t energy = digit->GetEnergy()/calibration;
 697     digit->SetEnergy(energy); //Now digit measures E in ADC counts
 698     Float_t time = digit->GetTime() ;
 699     time-=fcdb->GetTimeShiftEmc(module,column,row);
 700     digit->SetTime(time) ;
 701   }
 702 }
 703 //____________________________________________________________________________
 704 Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
 705   //Apply time calibration
 706   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 707
 708   //Determine rel.position of the cell absolute ID
 709   Int_t relId[4];
 710   geom->AbsToRelNumbering(absId,relId);
 711   Int_t module=relId[0];
 712   Int_t row   =relId[2];
 713   Int_t column=relId[3];
 714   time += fcdb->GetTimeShiftEmc(module,column,row);
 715   return time ;
 716 }
 717 //____________________________________________________________________________
 718 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
 719 {
 720   // Returns digitized value of the CPV charge in a pad absId
 721
 722   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 723
 724   //Determine rel.position of the cell absId
 725   Int_t relId[4];
 726   geom->AbsToRelNumbering(absId,relId);
 727   Int_t module=relId[0];
 728   Int_t row   =relId[2];
 729   Int_t column=relId[3];
 730
 731   Int_t channel = 0;
 732
 733   if(absId > fEmcCrystals){ //digitize CPV only
 734
 735     //reading calibration data for cell absId.
 736     Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
 737     Float_t adcChanelCpv   = fcdb->GetADCchannelCpv( module,column,row);
 738
 739     channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;
 740     Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
 741     if(channel > nMax ) channel = nMax ;
 742   }
 743   return channel ;
 744 }
 745
 746 //____________________________________________________________________________
 747 void AliPHOSDigitizer::Exec(Option_t *option)
 748 {
 749   // Steering method to process digitization for events
 750   // in the range from fFirstEvent to fLastEvent.
 751   // This range is optionally set by SetEventRange().
 752   // if fLastEvent=-1, then process events until the end.
 753   // by default fLastEvent = fFirstEvent (process only one event)
 754
 755   if (!fInit) { // to prevent overwrite existing file
 756     AliError(Form("Give a version name different from %s",
 757                   fEventFolderName.Data() )) ;
 758     return ;
 759   }
 760
 761   if (strstr(option,"print")) {
 762     Print();
 763     return ;
 764   }
 765
 766  if(strstr(option,"tim"))
 767      gBenchmark->Start("PHOSDigitizer");
 768
 769   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 770   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 771
 772   // Post Digitizer to the white board
 773   phosLoader->PostDigitizer(this) ;
 774
 775   if (fLastEvent == -1)
 776     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 777   else if (fManager)
 778     fLastEvent = fFirstEvent ;
 779
 780   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 781
 782   Int_t ievent ;
 783
 784   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 785     fEventCounter++ ;
 786
 787     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 788
 789     WriteDigits() ;
 790
 791     if(strstr(option,"deb"))
 792       PrintDigits(option);
 793
 794     //increment the total number of Digits per run
 795     fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;
 796  }
 797
 798   phosLoader->CleanDigitizer();
 799
 800   if(strstr(option,"tim")){
 801     gBenchmark->Stop("PHOSDigitizer");
 802     TString message ;
 803     message = "  took %f seconds for Digitizing %f seconds per event\n" ;
 804     AliInfo(Form( message.Data(),
 805          gBenchmark->GetCpuTime("PHOSDigitizer"),
 806          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
 807   }
 808 }
 809 //____________________________________________________________________________
 810 Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
 811   //calculate TOF resolution using beam-test resutls
 812   Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
 813   Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
 814   return TMath::Sqrt(a*a+b*b/e) ;
 815 }
 816
 817 ////____________________________________________________________________________
 818 //Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
 819 //{
 820 //  // Returns the shortest time among all time ticks
 821 //
 822 //  ticks->Sort() ; //Sort in accordance with times of ticks
 823 //  TIter it(ticks) ;
 824 //  AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
 825 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 826 //
 827 //  AliPHOSTick * t ;
 828 //  while((t=(AliPHOSTick*) it.Next())){
 829 //    if(t->GetTime() < time)  //This tick starts before crossing
 830 //      *ctick+=*t ;
 831 //    else
 832 //      return time ;
 833 //
 834 //    time = ctick->CrossingTime(fTimeThreshold) ;
 835 //  }
 836 //  return time ;
 837 //}
 838
 839 //____________________________________________________________________________
 840 Bool_t AliPHOSDigitizer::Init()
 841 {
 842   // Makes all memory allocations
 843   fInit = kTRUE ;
 844
 845   AliPHOSGeometry *geom;
 846   if (!(geom = AliPHOSGeometry::GetInstance()))
 847         geom = AliPHOSGeometry::GetInstance("IHEP","");
 848 //   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 849
 850   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
 851
 852   fFirstEvent = 0 ;
 853   fLastEvent = fFirstEvent ;
 854   if (fManager)
 855     fInput = fManager->GetNinputs() ;
 856   else
 857     fInput           = 1 ;
 858
 859   fInputFileNames  = new TString[fInput] ;
 860   fEventNames      = new TString[fInput] ;
 861   fInputFileNames[0] = GetTitle() ;
 862   fEventNames[0]     = fEventFolderName.Data() ;
 863   Int_t index ;
 864   for (index = 1 ; index < fInput ; index++) {
 865     fInputFileNames[index] = static_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 866     TString tempo = fManager->GetInputFolderName(index) ;
 867     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
 868   }
 869
 870   //to prevent cleaning of this object while GetEvent is called
 871   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 872   if(!rl){
 873     rl = AliRunLoader::Open(GetTitle(), fEventFolderName) ;
 874   }
 875   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 876   phosLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 877
 878   return fInit ;
 879 }
 880
 881 //____________________________________________________________________________
 882 void AliPHOSDigitizer::InitParameters()
 883 {
 884   // Set initial parameters Digitizer
 885
 886   fDigitsInRun  = 0 ;
 887   SetEventRange(0,-1) ;
 888   fPulse = new AliPHOSPulseGenerator();
 889   fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
 890   fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
 891
 892 }
 893
 894 //__________________________________________________________________
 895 void AliPHOSDigitizer::Print(const Option_t *)const
 896 {
 897   // Print Digitizer's parameters
 898   AliInfo(Form("\n------------------- %s -------------", GetName() )) ;
 899   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 900     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 901
 902     Int_t nStreams ;
 903     if (fManager)
 904       nStreams =  GetNInputStreams() ;
 905     else
 906       nStreams = fInput ;
 907
 908     Int_t index = 0 ;
 909     for (index = 0 ; index < nStreams ; index++) {
 910       TString tempo(fEventNames[index]) ;
 911       tempo += index ;
 912       printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ;
 913     }
 914
 915  //   printf("\nWith following parameters:\n") ;
 916  //   printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ;
 917  //   printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;
 918  //   printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ;
 919  //   printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;
 920     printf(" ---------------------------------------------------\n") ;
 921   }
 922   else
 923     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
 924
 925 }
 926
 927 //__________________________________________________________________
 928  void AliPHOSDigitizer::PrintDigits(Option_t * option)
 929 {
 930   // Print a table of digits
 931
 932
 933   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 934   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 935   TClonesArray * digits = phosLoader->Digits() ;
 936   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 937
 938   AliInfo(Form("%d", digits->GetEntriesFast())) ;
 939   printf("\nevent %d", gAlice->GetEvNumber()) ;
 940   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;
 941
 942
 943   if(strstr(option,"all")||strstr(option,"EMC")){
 944     //loop over digits
 945     AliPHOSDigit * digit;
 946     printf("\nEMC digits (with primaries):\n")  ;
 947     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 948     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
 949     Int_t index ;
 950     for (index = 0 ; (index < digits->GetEntriesFast()) &&
 951            (static_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
 952       digit = (AliPHOSDigit * )  digits->At(index) ;
 953       if(digit->GetNprimary() == 0)
 954         continue;
 955 //       printf("%6d  %8d    %6.5e %4d      %2d :",
 956 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
 957       printf("%6d  %.4f    %6.5e %4d      %2d :",
 958               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 959       Int_t iprimary;
 960       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 961         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 962       }
 963       printf("\n") ;
 964     }
 965   }
 966
 967   if(strstr(option,"all")||strstr(option,"CPV")){
 968
 969     //loop over CPV digits
 970     AliPHOSDigit * digit;
 971     printf("\nCPV digits (with primaries):\n")  ;
 972     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 973     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
 974     Int_t index ;
 975     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
 976       digit = (AliPHOSDigit * )  digits->At(index) ;
 977       if(digit->GetId() > maxEmc){
 978         printf("%6d  %8d    %4d      %2d :",
 979                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
 980         Int_t iprimary;
 981         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 982           printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 983         }
 984         printf("\n") ;
 985       }
 986     }
 987   }
 988
 989 }
 990
 991 //__________________________________________________________________
 992 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
 993 {  // Calculates the time signal generated by noise
 994   //PH  Info("TimeOfNoise", "Change me") ;
 995   return gRandom->Rndm() * 1.28E-5;
 996 }
 997
 998 //__________________________________________________________________
 999 void AliPHOSDigitizer::Unload()
1000 {
1001
1002   Int_t i ;
1003   for(i = 1 ; i < fInput ; i++){
1004     TString tempo(fEventNames[i]) ;
1005     tempo += i ;
1006     AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ;
1007     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1008     phosLoader->UnloadSDigits() ;
1009   }
1010
1011   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
1012   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1013   phosLoader->UnloadDigits() ;
1014 }
1015
1016 //____________________________________________________________________________
1017 void AliPHOSDigitizer::WriteDigits()
1018 {
1019
1020   // Makes TreeD in the output file.
1021   // Check if branch already exists:
1022   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1023   //      already existing branches.
1024   //   else creates branch with Digits, named "PHOS", title "...",
1025   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
1026   //      and names of files, from which digits are made.
1027
1028   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
1029   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1030
1031   const TClonesArray * digits = phosLoader->Digits() ;
1032   TTree * treeD = phosLoader->TreeD();
1033   if(!treeD){
1034     phosLoader->MakeTree("D");
1035     treeD = phosLoader->TreeD();
1036   }
1037
1038   // -- create Digits branch
1039   Int_t bufferSize = 32000 ;
1040   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
1041   digitsBranch->SetTitle(fEventFolderName);
1042   digitsBranch->Fill() ;
1043
1044   phosLoader->WriteDigits("OVERWRITE");
1045   phosLoader->WriteDigitizer("OVERWRITE");
1046
1047   Unload() ;
1048
1049 }