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(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     volpath = "/ALIC_1/PHOS_";
 287     volpath += i+1;
 288     if (gGeoManager->CheckPath(volpath.Data())) {
 289       isPresent[i]=1 ;
 290       nmod++ ;
 291     }
 292     else{
 293       isPresent[i]=0 ;
 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   //remove digits below thresholds
 584   Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
 585   for(Int_t i = 0 ; i < nEMC ; i++){
 586     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 587
 588     if(digit->GetEnergy() < emcThreshold){
 589       digits->RemoveAt(i) ;
 590       continue ;
 591     }
 592
 593     geom->AbsToRelNumbering(digit->GetId(),relId);
 594
 595 //    digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
 596
 597     Float_t tres = TimeResolution(digit->GetEnergy()) ;
 598     digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
 599
 600     fPulse->Reset();
 601     fPulse->SetAmplitude(digit->GetEnergy()/
 602                          fcdb->GetADCchannelEmc(relId[0],relId[3],relId[2]));
 603     fPulse->SetTZero(digit->GetTimeR());
 604     fPulse->MakeSamples();
 605     fPulse->GetSamples(fADCValuesHG, fADCValuesLG) ;
 606     Int_t nSamples = fPulse->GetRawFormatTimeBins();
 607     digit->SetALTROSamplesHG(nSamples,fADCValuesHG);
 608     digit->SetALTROSamplesLG(nSamples,fADCValuesLG);
 609   }
 610
 611   Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
 612   for(Int_t i = nEMC; i < nCPV ; i++){
 613     if( static_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
 614       digits->RemoveAt(i) ;
 615   }
 616
 617   digits->Compress() ;
 618   Int_t ndigits = digits->GetEntriesFast() ;
 619
 620   //Set indexes in list of digits and make true digitization of the energy
 621   for (Int_t i = 0 ; i < ndigits ; i++) {
 622     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 623     digit->SetIndexInList(i) ;
 624     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
 625       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
 626     }
 627   }
 628
 629 }
 630 //____________________________________________________________________________
 631 Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
 632   //Apply calibration
 633   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 634
 635   //Determine rel.position of the cell absolute ID
 636   Int_t relId[4];
 637   geom->AbsToRelNumbering(absId,relId);
 638   Int_t module=relId[0];
 639   Int_t row   =relId[2];
 640   Int_t column=relId[3];
 641   if(relId[1]==0){ //This Is EMC
 642     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
 643     return amp*calibration ;
 644   }
 645   return 0 ;
 646 }
 647 //____________________________________________________________________________
 648 void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
 649 {
 650   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
 651
 652   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 653
 654   //Determine rel.position of the cell absolute ID
 655   Int_t relId[4];
 656   geom->AbsToRelNumbering(digit->GetId(),relId);
 657   Int_t module=relId[0];
 658   Int_t row   =relId[2];
 659   Int_t column=relId[3];
 660   if(relId[1]==0){ //This Is EMC
 661     Float_t decalib     = fcdb->GetADCchannelEmcDecalib(module,column,row); // O(1)
 662     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row)*decalib;
 663     Float_t energy = digit->GetEnergy()/calibration;
 664     digit->SetEnergy(energy); //Now digit measures E in ADC counts
 665     Float_t time = digit->GetTime() ;
 666     time-=fcdb->GetTimeShiftEmc(module,column,row);
 667     digit->SetTime(time) ;
 668   }
 669 }
 670 //____________________________________________________________________________
 671 Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
 672   //Apply time calibration
 673   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 674
 675   //Determine rel.position of the cell absolute ID
 676   Int_t relId[4];
 677   geom->AbsToRelNumbering(absId,relId);
 678   Int_t module=relId[0];
 679   Int_t row   =relId[2];
 680   Int_t column=relId[3];
 681   time += fcdb->GetTimeShiftEmc(module,column,row);
 682   return time ;
 683 }
 684 //____________________________________________________________________________
 685 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
 686 {
 687   // Returns digitized value of the CPV charge in a pad absId
 688
 689   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 690
 691   //Determine rel.position of the cell absId
 692   Int_t relId[4];
 693   geom->AbsToRelNumbering(absId,relId);
 694   Int_t module=relId[0];
 695   Int_t row   =relId[2];
 696   Int_t column=relId[3];
 697
 698   Int_t channel = 0;
 699
 700   if(absId > fEmcCrystals){ //digitize CPV only
 701
 702     //reading calibration data for cell absId.
 703     Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
 704     Float_t adcChanelCpv   = fcdb->GetADCchannelCpv( module,column,row);
 705
 706     channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;
 707     Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
 708     if(channel > nMax ) channel = nMax ;
 709   }
 710   return channel ;
 711 }
 712
 713 //____________________________________________________________________________
 714 void AliPHOSDigitizer::Digitize(Option_t *option)
 715 {
 716   // Steering method to process digitization for events
 717   // in the range from fFirstEvent to fLastEvent.
 718   // This range is optionally set by SetEventRange().
 719   // if fLastEvent=-1, then process events until the end.
 720   // by default fLastEvent = fFirstEvent (process only one event)
 721
 722   if (!fInit) { // to prevent overwrite existing file
 723     AliError(Form("Give a version name different from %s",
 724                   fEventFolderName.Data() )) ;
 725     return ;
 726   }
 727
 728   if (strstr(option,"print")) {
 729     Print();
 730     return ;
 731   }
 732
 733  if(strstr(option,"tim"))
 734      gBenchmark->Start("PHOSDigitizer");
 735
 736   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 737   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 738
 739   if (fLastEvent == -1)
 740     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 741   else if (fDigInput)
 742     fLastEvent = fFirstEvent ;
 743
 744   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 745
 746   Int_t ievent ;
 747
 748   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 749     fEventCounter++ ;
 750
 751     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 752
 753     WriteDigits() ;
 754
 755     if(strstr(option,"deb"))
 756       PrintDigits(option);
 757
 758     //increment the total number of Digits per run
 759     fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;
 760  }
 761
 762   if(strstr(option,"tim")){
 763     gBenchmark->Stop("PHOSDigitizer");
 764     TString message ;
 765     message = "  took %f seconds for Digitizing %f seconds per event\n" ;
 766     AliInfo(Form( message.Data(),
 767          gBenchmark->GetCpuTime("PHOSDigitizer"),
 768          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
 769   }
 770 }
 771 //____________________________________________________________________________
 772 Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
 773   //calculate TOF resolution using beam-test resutls
 774   Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
 775   Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
 776   return TMath::Sqrt(a*a+b*b/e) ;
 777 }
 778
 779 ////____________________________________________________________________________
 780 //Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
 781 //{
 782 //  // Returns the shortest time among all time ticks
 783 //
 784 //  ticks->Sort() ; //Sort in accordance with times of ticks
 785 //  TIter it(ticks) ;
 786 //  AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
 787 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 788 //
 789 //  AliPHOSTick * t ;
 790 //  while((t=(AliPHOSTick*) it.Next())){
 791 //    if(t->GetTime() < time)  //This tick starts before crossing
 792 //      *ctick+=*t ;
 793 //    else
 794 //      return time ;
 795 //
 796 //    time = ctick->CrossingTime(fTimeThreshold) ;
 797 //  }
 798 //  return time ;
 799 //}
 800
 801 //____________________________________________________________________________
 802 Bool_t AliPHOSDigitizer::Init()
 803 {
 804   // Makes all memory allocations
 805   fInit = kTRUE ;
 806
 807   AliPHOSGeometry *geom;
 808   if (!(geom = AliPHOSGeometry::GetInstance()))
 809         geom = AliPHOSGeometry::GetInstance("IHEP","");
 810 //   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 811
 812   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
 813
 814   fFirstEvent = 0 ;
 815   fLastEvent = fFirstEvent ;
 816   if (fDigInput)
 817     fInput = fDigInput->GetNinputs() ;
 818   else
 819     fInput           = 1 ;
 820
 821   fInputFileNames  = new TString[fInput] ;
 822   fEventNames      = new TString[fInput] ;
 823   fInputFileNames[0] = GetTitle() ;
 824   fEventNames[0]     = fEventFolderName.Data() ;
 825   Int_t index ;
 826   for (index = 1 ; index < fInput ; index++) {
 827     fInputFileNames[index] = static_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
 828     TString tempo = fDigInput->GetInputFolderName(index) ;
 829     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fDigInput
 830   }
 831
 832   //to prevent cleaning of this object while GetEvent is called
 833   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 834   if(!rl){
 835     AliRunLoader::Open(GetTitle(), fEventFolderName) ;
 836   }
 837   return fInit ;
 838 }
 839
 840 //____________________________________________________________________________
 841 void AliPHOSDigitizer::InitParameters()
 842 {
 843   // Set initial parameters Digitizer
 844
 845   fDigitsInRun  = 0 ;
 846   SetEventRange(0,-1) ;
 847   fPulse = new AliPHOSPulseGenerator();
 848   fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
 849   fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
 850
 851 }
 852
 853 //__________________________________________________________________
 854 void AliPHOSDigitizer::Print(const Option_t *)const
 855 {
 856   // Print Digitizer's parameters
 857   AliInfo(Form("\n------------------- %s -------------", GetName() )) ;
 858   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 859     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 860
 861     Int_t nStreams ;
 862     if (fDigInput)
 863       nStreams =  GetNInputStreams() ;
 864     else
 865       nStreams = fInput ;
 866
 867     Int_t index = 0 ;
 868     for (index = 0 ; index < nStreams ; index++) {
 869       TString tempo(fEventNames[index]) ;
 870       tempo += index ;
 871       printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ;
 872     }
 873
 874  //   printf("\nWith following parameters:\n") ;
 875  //   printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ;
 876  //   printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;
 877  //   printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ;
 878  //   printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;
 879     printf(" ---------------------------------------------------\n") ;
 880   }
 881   else
 882     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
 883
 884 }
 885
 886 //__________________________________________________________________
 887  void AliPHOSDigitizer::PrintDigits(Option_t * option)
 888 {
 889   // Print a table of digits
 890
 891
 892   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 893   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 894   TClonesArray * digits = phosLoader->Digits() ;
 895   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 896
 897   AliInfo(Form("%d", digits->GetEntriesFast())) ;
 898   printf("\nevent %d", gAlice->GetEvNumber()) ;
 899   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;
 900
 901
 902   if(strstr(option,"all")||strstr(option,"EMC")){
 903     //loop over digits
 904     AliPHOSDigit * digit;
 905     printf("\nEMC digits (with primaries):\n")  ;
 906     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 907     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
 908     Int_t index ;
 909     for (index = 0 ; (index < digits->GetEntriesFast()) &&
 910            (static_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
 911       digit = (AliPHOSDigit * )  digits->At(index) ;
 912       if(digit->GetNprimary() == 0)
 913         continue;
 914 //       printf("%6d  %8d    %6.5e %4d      %2d :",
 915 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
 916       printf("%6d  %.4f    %6.5e %4d      %2d :",
 917               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 918       Int_t iprimary;
 919       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 920         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 921       }
 922       printf("\n") ;
 923     }
 924   }
 925
 926   if(strstr(option,"all")||strstr(option,"CPV")){
 927
 928     //loop over CPV digits
 929     AliPHOSDigit * digit;
 930     printf("\nCPV digits (with primaries):\n")  ;
 931     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 932     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
 933     Int_t index ;
 934     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
 935       digit = (AliPHOSDigit * )  digits->At(index) ;
 936       if(digit->GetId() > maxEmc){
 937         printf("%6d  %8d    %4d      %2d :",
 938                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
 939         Int_t iprimary;
 940         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 941           printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 942         }
 943         printf("\n") ;
 944       }
 945     }
 946   }
 947
 948 }
 949
 950 //__________________________________________________________________
 951 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
 952 {  // Calculates the time signal generated by noise
 953   //PH  Info("TimeOfNoise", "Change me") ;
 954   return gRandom->Rndm() * 1.28E-5;
 955 }
 956
 957 //__________________________________________________________________
 958 void AliPHOSDigitizer::Unload()
 959 {
 960
 961   Int_t i ;
 962   for(i = 1 ; i < fInput ; i++){
 963     TString tempo(fEventNames[i]) ;
 964     tempo += i ;
 965     AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ;
 966     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 967     phosLoader->UnloadSDigits() ;
 968   }
 969
 970   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 971   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 972   phosLoader->UnloadDigits() ;
 973 }
 974
 975 //____________________________________________________________________________
 976 void AliPHOSDigitizer::WriteDigits()
 977 {
 978
 979   // Makes TreeD in the output file.
 980   // Check if branch already exists:
 981   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
 982   //      already existing branches.
 983   //   else creates branch with Digits, named "PHOS", title "...",
 984   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
 985   //      and names of files, from which digits are made.
 986
 987   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 988   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 989
 990   const TClonesArray * digits = phosLoader->Digits() ;
 991   TTree * treeD = phosLoader->TreeD();
 992   if(!treeD){
 993     phosLoader->MakeTree("D");
 994     treeD = phosLoader->TreeD();
 995   }
 996
 997   // -- create Digits branch
 998   Int_t bufferSize = 32000 ;
 999   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
1000   digitsBranch->SetTitle(fEventFolderName);
1001   digitsBranch->Fill() ;
1002
1003   phosLoader->WriteDigits("OVERWRITE");
1004
1005   Unload() ;
1006
1007 }