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