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