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