PHOS/AliPHOSDigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /* History of cvs commits:
  19  *
  20  * $Log: AliPHOSDigitizer.cxx,v $
  21  * Revision 1.104  2007/12/18 09:08:18  hristov
  22  * Splitting of the QA maker into simulation and reconstruction dependent parts (Yves)
  23  *
  24  * Revision 1.103  2007/11/07 11:25:06  schutz
  25  * Comment out the QA checking before starting digitization
  26  *
  27  * Revision 1.102  2007/10/19 18:04:29  schutz
  28  * The standalone QA data maker is called from AliSimulation and AliReconstruction outside the event loop; i.e. re-reading the data. The QA data making in the event loop has been commented out.
  29  *
  30  * Revision 1.101  2007/10/14 21:08:10  schutz
  31  * Introduced the checking of QA results from previous step before entering the event loop
  32  *
  33  * Revision 1.100  2007/10/10 09:05:10  schutz
  34  * Changing name QualAss to QA
  35  *
  36  * Revision 1.99  2007/09/30 17:08:20  schutz
  37  * Introducing the notion of QA data acquisition cycle (needed by online)
  38  *
  39  * Revision 1.98  2007/09/26 14:22:17  cvetan
  40  * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
  41  *
  42  * Revision 1.97  2007/08/07 14:12:03  kharlov
  43  * Quality assurance added (Yves Schutz)
  44  *
  45  * Revision 1.96  2007/04/28 10:43:36  policheh
  46  * Dead channels simulation: digit energy sets to 0.
  47  *
  48  * Revision 1.95  2007/04/10 07:20:52  kharlov
  49  * Decalibration should use the same CDB as calibration in AliPHOSClusterizerv1
  50  *
  51  * Revision 1.94  2007/02/01 10:34:47  hristov
  52  * Removing warnings on Solaris x86
  53  *
  54  * Revision 1.93  2006/10/17 13:17:01  kharlov
  55  * Replace AliInfo by AliDebug
  56  *
  57  * Revision 1.92  2006/08/28 10:01:56  kharlov
  58  * Effective C++ warnings fixed (Timur Pocheptsov)
  59  *
  60  * Revision 1.91  2006/04/29 20:25:30  hristov
  61  * Decalibration is implemented (Yu.Kharlov)
  62  *
  63  * Revision 1.90  2006/04/22 10:30:17  hristov
  64  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
  65  *
  66  * Revision 1.89  2006/04/11 15:22:59  hristov
  67  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
  68  *
  69  * Revision 1.88  2006/03/13 14:05:43  kharlov
  70  * Calibration objects for EMC and CPV
  71  *
  72  * Revision 1.87  2005/08/24 15:33:49  kharlov
  73  * Calibration data for raw digits
  74  *
  75  * Revision 1.86  2005/07/12 20:07:35  hristov
  76  * Changes needed to run simulation and reconstrruction in the same AliRoot session
  77  *
  78  * Revision 1.85  2005/05/28 14:19:04  schutz
  79  * Compilation warnings fixed by T.P.
  80  *
  81  */
  82
  83 //_________________________________________________________________________
  84 //*-- Author :  Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  85 //////////////////////////////////////////////////////////////////////////////
  86 // This TTask performs digitization of Summable digits (in the PHOS case it is just
  87 // the sum of contributions from all primary particles into a given cell).
  88 // In addition it performs mixing of summable digits from different events.
  89 // The name of the TTask is also the title of the branch that will contain
  90 // the created SDigits
  91 // The title of the TTAsk is the name of the file that contains the hits from
  92 // which the SDigits are created
  93 //
  94 // For each event two branches are created in TreeD:
  95 //   "PHOS" - list of digits
  96 //   "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
  97 //
  98 // Note, that one can set a title for new digits branch, and repeat digitization with
  99 // another set of parameters.
 100 //
 101 // Use case:
 102 // root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
 103 // root[1] d->ExecuteTask()
 104 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
 105 //                       //Digitizes SDigitis in all events found in file galice.root
 106 //
 107 // root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;
 108 //                       // Will read sdigits from galice1.root
 109 // root[3] d1->MixWith("galice2.root")
 110 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
 111 //                       // Reads another set of sdigits from galice2.root
 112 // root[3] d1->MixWith("galice3.root")
 113 //                       // Reads another set of sdigits from galice3.root
 114 // root[4] d->ExecuteTask("deb timing")
 115 //                       // Reads SDigits from files galice1.root, galice2.root ....
 116 //                       // mixes them and stores produced Digits in file galice1.root
 117 //                       // deb - prints number of produced digits
 118 //                       // deb all - prints list of produced digits
 119 //                       // timing  - prints time used for digitization
 120 //
 121
 122 // --- ROOT system ---
 123 #include "TTree.h"
 124 #include "TSystem.h"
 125 #include "TBenchmark.h"
 126 #include "TRandom.h"
 127 #include "TMath.h"
 128
 129 // --- Standard library ---
 130
 131 // --- AliRoot header files ---
 132 #include <TGeoManager.h>
 133 #include "AliLog.h"
 134 #include "AliRunDigitizer.h"
 135 #include "AliPHOSDigit.h"
 136 #include "AliPHOSDigitizer.h"
 137 #include "AliPHOSGeometry.h"
 138 #include "AliPHOSTick.h"
 139 #include "AliPHOSSimParam.h"
 140 #include "AliPHOSCalibData.h"
 141 #include "AliRunLoader.h"
 142 #include "AliPHOSLoader.h"
 143 #include "AliPHOSPulseGenerator.h"
 144
 145 ClassImp(AliPHOSDigitizer)
 146
 147
 148 //____________________________________________________________________________
 149 AliPHOSDigitizer::AliPHOSDigitizer() :
 150   AliDigitizer("",""),
 151   fDefaultInit(kTRUE),
 152   fDigitsInRun(0),
 153   fInit(kFALSE),
 154   fInput(0),
 155   fInputFileNames(0x0),
 156   fEventNames(0x0),
 157   fEmcCrystals(0),
 158   fEventFolderName(""),
 159   fFirstEvent(0),
 160   fLastEvent(0),
 161   fcdb(0x0),
 162   fEventCounter(0),
 163   fPulse(0),
 164   fADCValuesLG(0),
 165   fADCValuesHG(0)
 166 {
 167   // ctor
 168   InitParameters() ;
 169   fManager = 0 ;                     // We work in the standalong mode
 170 }
 171
 172 //____________________________________________________________________________
 173 AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName,
 174                                    TString eventFolderName):
 175   AliDigitizer("PHOS"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
 176   fDefaultInit(kFALSE),
 177   fDigitsInRun(0),
 178   fInit(kFALSE),
 179   fInput(0),
 180   fInputFileNames(0x0),
 181   fEventNames(0x0),
 182   fEmcCrystals(0),
 183   fEventFolderName(eventFolderName),
 184   fFirstEvent(0),
 185   fLastEvent(0),
 186   fcdb(0x0),
 187   fEventCounter(0),
 188   fPulse(0),
 189   fADCValuesLG(0),
 190   fADCValuesHG(0)
 191 {
 192   // ctor
 193   InitParameters() ;
 194   Init() ;
 195   fDefaultInit = kFALSE ;
 196   fManager = 0 ;                     // We work in the standalone mode
 197   fcdb = new AliPHOSCalibData(-1);
 198 }
 199
 200 //____________________________________________________________________________
 201 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) :
 202   AliDigitizer(d),
 203   fDefaultInit(d.fDefaultInit),
 204   fDigitsInRun(d.fDigitsInRun),
 205   fInit(d.fInit),
 206   fInput(d.fInput),
 207   fInputFileNames(0x0),//?
 208   fEventNames(0x0),//?
 209   fEmcCrystals(d.fEmcCrystals),
 210   fEventFolderName(d.fEventFolderName),
 211   fFirstEvent(d.fFirstEvent),
 212   fLastEvent(d.fLastEvent),
 213   fcdb (0x0),
 214   fEventCounter(0),
 215   fPulse(0),
 216   fADCValuesLG(0),
 217   fADCValuesHG(0)
 218 {
 219   // copyy ctor
 220   SetName(d.GetName()) ;
 221   SetTitle(d.GetTitle()) ;
 222   fcdb = new AliPHOSCalibData(-1);
 223 }
 224
 225 //____________________________________________________________________________
 226 AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd) :
 227   AliDigitizer(rd,"PHOS"+AliConfig::Instance()->GetDigitizerTaskName()),
 228   fDefaultInit(kFALSE),
 229   fDigitsInRun(0),
 230   fInit(kFALSE),
 231   fInput(0),
 232   fInputFileNames(0x0),
 233   fEventNames(0x0),
 234   fEmcCrystals(0),
 235   fEventFolderName(fManager->GetInputFolderName(0)),
 236   fFirstEvent(0),
 237   fLastEvent(0),
 238   fcdb (0x0),
 239   fEventCounter(0),
 240   fPulse(0),
 241   fADCValuesLG(0),
 242   fADCValuesHG(0)
 243
 244 {
 245   // ctor Init() is called by RunDigitizer
 246   fManager = rd ;
 247   SetTitle(static_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
 248   InitParameters() ;
 249   fDefaultInit = kFALSE ;
 250   fcdb = new AliPHOSCalibData(-1);
 251 }
 252
 253 //____________________________________________________________________________
 254   AliPHOSDigitizer::~AliPHOSDigitizer()
 255 {
 256   // dtor
 257   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 258   if(rl){
 259     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 260
 261     if(phosLoader)
 262       phosLoader->CleanDigitizer() ;
 263   }
 264
 265   delete [] fInputFileNames ;
 266   delete [] fEventNames ;
 267
 268   delete fPulse;
 269   delete [] fADCValuesLG;
 270   delete [] fADCValuesHG;
 271
 272   if(fcdb){ delete fcdb ; fcdb=0;}
 273
 274 }
 275
 276 //____________________________________________________________________________
 277 void AliPHOSDigitizer::Digitize(Int_t event)
 278 {
 279
 280   // Makes the digitization of the collected summable digits.
 281   //  It first creates the array of all PHOS modules
 282   //  filled with noise (different for EMC, and CPV) and
 283   //  then adds contributions from SDigits.
 284   // This design avoids scanning over the list of digits to add
 285   // contribution to new SDigits only.
 286
 287   Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
 288
 289   //First stream
 290   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 291   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 292
 293   Int_t readEvent = event ;
 294   if (fManager)
 295     readEvent = static_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
 296   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 297                   readEvent, GetTitle(), fEventFolderName.Data())) ;
 298   rl->GetEvent(readEvent) ;
 299   phosLoader->CleanSDigits() ;
 300   phosLoader->LoadSDigits("READ") ;
 301
 302   //Prepare Output
 303   TClonesArray * digits = phosLoader->Digits() ;
 304   if( !digits ) {
 305     phosLoader->MakeDigitsArray() ;
 306     digits = phosLoader->Digits() ;
 307   }
 308
 309   digits->Clear() ;
 310
 311   //
 312   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 313   //Making digits with noise, first EMC
 314   //Check which PHOS modules are present
 315   Bool_t isPresent[5] ;
 316   TString volpath ;
 317   Int_t nmod=0 ;
 318   for(Int_t i=0; i<5; i++){
 319     volpath = "/ALIC_1/PHOS_";
 320     volpath += i+1;
 321     if (gGeoManager->CheckPath(volpath.Data())) {
 322       isPresent[i]=1 ;
 323       nmod++ ;
 324     }
 325     else{
 326       isPresent[i]=0 ;
 327     }
 328   }
 329
 330   Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
 331
 332   Int_t nCPV ;
 333   Int_t absID ;
 334
 335   //check if CPV exists
 336   Bool_t isCPVpresent=0 ;
 337   for(Int_t i=1; i<=5 && !isCPVpresent; i++){
 338     volpath = "/ALIC_1/PHOS_";
 339     volpath += i;
 340     volpath += "/PCPV_1";
 341     if (gGeoManager->CheckPath(volpath.Data()))
 342       isCPVpresent=1 ;
 343   }
 344
 345   if(isCPVpresent){
 346     nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
 347   }
 348   else{
 349      nCPV = nEMC ;
 350   }
 351
 352   digits->Expand(nCPV) ;
 353
 354   //take all the inputs to add together and load the SDigits
 355   TObjArray * sdigArray = new TObjArray(fInput) ;
 356   sdigArray->AddAt(phosLoader->SDigits(), 0) ;
 357
 358   for(Int_t i = 1 ; i < fInput ; i++){
 359     TString tempo(fEventNames[i]) ;
 360     tempo += i ;
 361     AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
 362     if(!rl2){
 363       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
 364       if (!rl2) {
 365         Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
 366         return ;
 367       }
 368       rl2->LoadHeader();
 369     }
 370     AliPHOSLoader * phosLoader2 = static_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
 371
 372     if(fManager){
 373       readEvent = static_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
 374     }
 375     TClonesArray * digs ;
 376     if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
 377       AliInfo(Form("Adding event %d from input stream %d %s %s",
 378                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
 379        rl2->GetEvent(readEvent) ;
 380        phosLoader2->CleanDigits() ;
 381        phosLoader2->LoadDigits("READ") ;
 382        digs = phosLoader2->Digits() ;
 383        toMakeNoise=0 ; //Do not add noise, it is already in stream
 384     }
 385     else{
 386       AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
 387                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
 388        rl2->GetEvent(readEvent) ;
 389        phosLoader2->CleanSDigits() ;
 390        phosLoader2->LoadSDigits("READ") ;
 391        digs = phosLoader2->SDigits() ;
 392     }
 393     sdigArray->AddAt(digs, i) ;
 394   }
 395
 396   //Find the first crystal with signal
 397   Int_t nextSig = 200000 ;
 398   TClonesArray * sdigits ;
 399   for(Int_t i = 0 ; i < fInput ; i++){
 400     sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 401     if ( !sdigits->GetEntriesFast() )
 402       continue ;
 403     Int_t curNext = static_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
 404     if(curNext < nextSig)
 405       nextSig = curNext ;
 406   }
 407
 408   TArrayI index(fInput) ;
 409   index.Reset() ;  //Set all indexes to zero
 410
 411   AliPHOSDigit * digit ;
 412   AliPHOSDigit * curSDigit ;
 413
 414 //  TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
 415
 416   //Put Noise contribution
 417   Double_t apdNoise = 0. ;
 418   if(toMakeNoise)
 419      apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ;
 420
 421   Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
 422   Int_t idigit= 0;
 423   for(Int_t imod=0; imod<5; imod++){
 424     if(!isPresent[imod])
 425       continue ;
 426     Int_t firstAbsId=imod*emcpermod+1 ;
 427     Int_t lastAbsId =(imod+1)*emcpermod ;
 428     for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
 429       Float_t noise = gRandom->Gaus(0.,apdNoise) ;
 430       new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
 431       //look if we have to add signal?
 432       digit = static_cast<AliPHOSDigit *>(digits->At(idigit)) ;
 433       idigit++ ;
 434
 435       if(absID==nextSig){
 436       //Add SDigits from all inputs
 437 //      ticks->Clear() ;
 438 //      Int_t contrib = 0 ;
 439
 440 //New Timing model is necessary
 441 //      Float_t a = digit->GetEnergy() ;
 442 //      Float_t b = TMath::Abs( a / fTimeSignalLength) ;
 443 //      //Mark the beginning of the signal
 444 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);
 445 //      //Mark the end of the signal
 446 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b);
 447
 448 // Calculate time as time of the largest digit
 449         Float_t time = digit->GetTime() ;
 450         Float_t eTime= digit->GetEnergy() ;
 451
 452         //loop over inputs
 453         for(Int_t i = 0 ; i < fInput ; i++){
 454           if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
 455             curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 456             if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
 457               curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
 458               curSDigit->SetTime(CalibrateT(curSDigit->GetTime(),curSDigit->GetId())) ;
 459             }
 460           }
 461           else
 462             curSDigit = 0 ;
 463           //May be several digits will contribute from the same input
 464           while(curSDigit && curSDigit->GetId() == absID){
 465             //Shift primary to separate primaries belonging different inputs
 466             Int_t primaryoffset ;
 467             if(fManager)
 468               primaryoffset = fManager->GetMask(i) ;
 469             else
 470               primaryoffset = 10000000*i ;
 471             curSDigit->ShiftPrimary(primaryoffset) ;
 472
 473 //New Timing model is necessary
 474 //        a = curSDigit->GetEnergy() ;
 475 //        b = a /fTimeSignalLength ;
 476 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);
 477 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
 478           if(curSDigit->GetEnergy()>eTime){
 479             eTime=curSDigit->GetEnergy() ;
 480             time=curSDigit->GetTime() ;
 481           }
 482             *digit += *curSDigit ;  //add energies
 483
 484             index[i]++ ;
 485             if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 486               curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 487             else
 488               curSDigit = 0 ;
 489           }
 490         }
 491
 492         //calculate and set time
 493 //New Timing model is necessary
 494 //      Float_t time = FrontEdgeTime(ticks) ;
 495         digit->SetTime(time) ;
 496
 497         //Find next signal module
 498         nextSig = 200000 ;
 499         for(Int_t i = 0 ; i < fInput ; i++){
 500           sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 501           Int_t curNext = nextSig ;
 502           if(sdigits->GetEntriesFast() > index[i] ){
 503             curNext = static_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
 504           }
 505           if(curNext < nextSig) nextSig = curNext ;
 506         }
 507       }
 508     }
 509   }
 510
 511   //distretize energy if necessary
 512   if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
 513     Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
 514     for(Int_t i = 0 ; i < nEMC ; i++){
 515       digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 516       digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
 517     }
 518   }
 519   //Apply decalibration if necessary
 520   for(Int_t i = 0 ; i < nEMC ; i++){
 521     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 522     Decalibrate(digit) ;
 523   }
 524
 525
 526 //  ticks->Delete() ;
 527 //  delete ticks ;
 528
 529   //Now CPV digits (different noise and no timing)
 530   Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
 531   Int_t nEMCtotal=emcpermod*5 ;
 532   Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
 533   if(isCPVpresent){  //CPV is present in current geometry
 534     for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
 535       if(!isPresent[imod])
 536         continue ;
 537       Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
 538       Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
 539       for(absID = firstAbsId; absID <= lastAbsId; absID++){
 540         Float_t noise = gRandom->Gaus(0., cpvNoise) ;
 541         new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
 542         idigit++ ;
 543         //look if we have to add signal?
 544         if(absID==nextSig){
 545           digit = static_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
 546           //Add SDigits from all inputs
 547           for(Int_t i = 0 ; i < fInput ; 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              //May be several digits will contribute from the same input
 554              while(curSDigit && curSDigit->GetId() == absID){
 555                //Shift primary to separate primaries belonging different inputs
 556                Int_t primaryoffset ;
 557                if(fManager)
 558                  primaryoffset = fManager->GetMask(i) ;
 559                else
 560                  primaryoffset = 10000000*i ;
 561                curSDigit->ShiftPrimary(primaryoffset) ;
 562
 563                //add energies
 564                *digit += *curSDigit ;
 565                index[i]++ ;
 566                if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 567                  curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;
 568                else
 569                  curSDigit = 0 ;
 570              }
 571           }
 572
 573           //Find next signal module
 574           nextSig = 200000 ;
 575           for(Int_t i = 0 ; i < fInput ; i++){
 576             sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 577             Int_t curNext = nextSig ;
 578             if(sdigits->GetEntriesFast() > index[i] )
 579               curNext = static_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
 580             if(curNext < nextSig) nextSig = curNext ;
 581           }
 582
 583         }
 584       }
 585     }
 586   }
 587
 588   delete sdigArray ; //We should not delete its contents
 589
 590   Int_t relId[4];
 591
 592   //set amplitudes in bad channels to zero
 593
 594   for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
 595     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 596     geom->AbsToRelNumbering(digit->GetId(),relId);
 597     if(relId[1] == 0) // Emc
 598       if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.);
 599   }
 600
 601   //remove digits below thresholds
 602   Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
 603   for(Int_t i = 0 ; i < nEMC ; i++){
 604     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 605
 606     if(digit->GetEnergy() < emcThreshold){
 607       digits->RemoveAt(i) ;
 608       continue ;
 609     }
 610
 611     geom->AbsToRelNumbering(digit->GetId(),relId);
 612
 613     digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
 614
 615     Float_t tres = TimeResolution(digit->GetEnergy()) ;
 616     digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
 617
 618     fPulse->Reset();
 619     fPulse->SetAmplitude(digit->GetEnergy()/
 620                          fcdb->GetADCchannelEmc(relId[0],relId[3],relId[2]));
 621     fPulse->SetTZero(digit->GetTimeR());
 622     fPulse->MakeSamples();
 623     fPulse->GetSamples(fADCValuesHG, fADCValuesLG) ;
 624     Int_t nSamples = fPulse->GetRawFormatTimeBins();
 625     digit->SetALTROSamplesHG(nSamples,fADCValuesHG);
 626     digit->SetALTROSamplesLG(nSamples,fADCValuesLG);
 627   }
 628
 629   Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
 630   for(Int_t i = nEMC; i < nCPV ; i++){
 631     if( static_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
 632       digits->RemoveAt(i) ;
 633   }
 634
 635   digits->Compress() ;
 636   Int_t ndigits = digits->GetEntriesFast() ;
 637
 638   //Set indexes in list of digits and make true digitization of the energy
 639   for (Int_t i = 0 ; i < ndigits ; i++) {
 640     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 641     digit->SetIndexInList(i) ;
 642     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
 643       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
 644     }
 645   }
 646
 647 }
 648 //____________________________________________________________________________
 649 Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
 650   //Apply calibration
 651   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 652
 653   //Determine rel.position of the cell absolute ID
 654   Int_t relId[4];
 655   geom->AbsToRelNumbering(absId,relId);
 656   Int_t module=relId[0];
 657   Int_t row   =relId[2];
 658   Int_t column=relId[3];
 659   if(relId[1]==0){ //This Is EMC
 660     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
 661     return amp*calibration ;
 662   }
 663   return 0 ;
 664 }
 665 //____________________________________________________________________________
 666 void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
 667 {
 668   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
 669
 670   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 671
 672   //Determine rel.position of the cell absolute ID
 673   Int_t relId[4];
 674   geom->AbsToRelNumbering(digit->GetId(),relId);
 675   Int_t module=relId[0];
 676   Int_t row   =relId[2];
 677   Int_t column=relId[3];
 678   if(relId[1]==0){ //This Is EMC
 679     Float_t decalib     = fcdb->GetADCchannelEmcDecalib(module,column,row); // O(1)
 680     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row)*decalib;
 681     Float_t energy = digit->GetEnergy()/calibration;
 682     digit->SetEnergy(energy); //Now digit measures E in ADC counts
 683     Float_t time = digit->GetTime() ;
 684     time-=fcdb->GetTimeShiftEmc(module,column,row);
 685     digit->SetTime(time) ;
 686   }
 687 }
 688 //____________________________________________________________________________
 689 Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
 690   //Apply time calibration
 691   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 692
 693   //Determine rel.position of the cell absolute ID
 694   Int_t relId[4];
 695   geom->AbsToRelNumbering(absId,relId);
 696   Int_t module=relId[0];
 697   Int_t row   =relId[2];
 698   Int_t column=relId[3];
 699   time += fcdb->GetTimeShiftEmc(module,column,row);
 700   return time ;
 701 }
 702 //____________________________________________________________________________
 703 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
 704 {
 705   // Returns digitized value of the CPV charge in a pad absId
 706
 707   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 708
 709   //Determine rel.position of the cell absId
 710   Int_t relId[4];
 711   geom->AbsToRelNumbering(absId,relId);
 712   Int_t module=relId[0];
 713   Int_t row   =relId[2];
 714   Int_t column=relId[3];
 715
 716   Int_t channel = 0;
 717
 718   if(absId > fEmcCrystals){ //digitize CPV only
 719
 720     //reading calibration data for cell absId.
 721     Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
 722     Float_t adcChanelCpv   = fcdb->GetADCchannelCpv( module,column,row);
 723
 724     channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;
 725     Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
 726     if(channel > nMax ) channel = nMax ;
 727   }
 728   return channel ;
 729 }
 730
 731 //____________________________________________________________________________
 732 void AliPHOSDigitizer::Exec(Option_t *option)
 733 {
 734   // Steering method to process digitization for events
 735   // in the range from fFirstEvent to fLastEvent.
 736   // This range is optionally set by SetEventRange().
 737   // if fLastEvent=-1, then process events until the end.
 738   // by default fLastEvent = fFirstEvent (process only one event)
 739
 740   if (!fInit) { // to prevent overwrite existing file
 741     AliError(Form("Give a version name different from %s",
 742                   fEventFolderName.Data() )) ;
 743     return ;
 744   }
 745
 746   if (strstr(option,"print")) {
 747     Print();
 748     return ;
 749   }
 750
 751  if(strstr(option,"tim"))
 752      gBenchmark->Start("PHOSDigitizer");
 753
 754   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 755   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 756
 757   // Post Digitizer to the white board
 758   phosLoader->PostDigitizer(this) ;
 759
 760   if (fLastEvent == -1)
 761     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 762   else if (fManager)
 763     fLastEvent = fFirstEvent ;
 764
 765   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 766
 767   Int_t ievent ;
 768
 769   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 770     fEventCounter++ ;
 771
 772     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 773
 774     WriteDigits() ;
 775
 776     if(strstr(option,"deb"))
 777       PrintDigits(option);
 778
 779     //increment the total number of Digits per run
 780     fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;
 781  }
 782
 783   phosLoader->CleanDigitizer();
 784
 785   if(strstr(option,"tim")){
 786     gBenchmark->Stop("PHOSDigitizer");
 787     TString message ;
 788     message = "  took %f seconds for Digitizing %f seconds per event\n" ;
 789     AliInfo(Form( message.Data(),
 790          gBenchmark->GetCpuTime("PHOSDigitizer"),
 791          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
 792   }
 793 }
 794 //____________________________________________________________________________
 795 Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
 796   //calculate TOF resolution using beam-test resutls
 797   Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
 798   Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
 799   return TMath::Sqrt(a*a+b*b/e) ;
 800 }
 801
 802 ////____________________________________________________________________________
 803 //Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
 804 //{
 805 //  // Returns the shortest time among all time ticks
 806 //
 807 //  ticks->Sort() ; //Sort in accordance with times of ticks
 808 //  TIter it(ticks) ;
 809 //  AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
 810 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 811 //
 812 //  AliPHOSTick * t ;
 813 //  while((t=(AliPHOSTick*) it.Next())){
 814 //    if(t->GetTime() < time)  //This tick starts before crossing
 815 //      *ctick+=*t ;
 816 //    else
 817 //      return time ;
 818 //
 819 //    time = ctick->CrossingTime(fTimeThreshold) ;
 820 //  }
 821 //  return time ;
 822 //}
 823
 824 //____________________________________________________________________________
 825 Bool_t AliPHOSDigitizer::Init()
 826 {
 827   // Makes all memory allocations
 828   fInit = kTRUE ;
 829
 830   AliPHOSGeometry *geom;
 831   if (!(geom = AliPHOSGeometry::GetInstance()))
 832         geom = AliPHOSGeometry::GetInstance("IHEP","");
 833 //   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 834
 835   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
 836
 837   fFirstEvent = 0 ;
 838   fLastEvent = fFirstEvent ;
 839   if (fManager)
 840     fInput = fManager->GetNinputs() ;
 841   else
 842     fInput           = 1 ;
 843
 844   fInputFileNames  = new TString[fInput] ;
 845   fEventNames      = new TString[fInput] ;
 846   fInputFileNames[0] = GetTitle() ;
 847   fEventNames[0]     = fEventFolderName.Data() ;
 848   Int_t index ;
 849   for (index = 1 ; index < fInput ; index++) {
 850     fInputFileNames[index] = static_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 851     TString tempo = fManager->GetInputFolderName(index) ;
 852     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
 853   }
 854
 855   //to prevent cleaning of this object while GetEvent is called
 856   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 857   if(!rl){
 858     rl = AliRunLoader::Open(GetTitle(), fEventFolderName) ;
 859   }
 860   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 861   phosLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 862
 863   return fInit ;
 864 }
 865
 866 //____________________________________________________________________________
 867 void AliPHOSDigitizer::InitParameters()
 868 {
 869   // Set initial parameters Digitizer
 870
 871   fDigitsInRun  = 0 ;
 872   SetEventRange(0,-1) ;
 873   fPulse = new AliPHOSPulseGenerator();
 874   fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
 875   fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
 876
 877 }
 878
 879 //__________________________________________________________________
 880 void AliPHOSDigitizer::Print(const Option_t *)const
 881 {
 882   // Print Digitizer's parameters
 883   AliInfo(Form("\n------------------- %s -------------", GetName() )) ;
 884   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 885     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 886
 887     Int_t nStreams ;
 888     if (fManager)
 889       nStreams =  GetNInputStreams() ;
 890     else
 891       nStreams = fInput ;
 892
 893     Int_t index = 0 ;
 894     for (index = 0 ; index < nStreams ; index++) {
 895       TString tempo(fEventNames[index]) ;
 896       tempo += index ;
 897       printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ;
 898     }
 899
 900  //   printf("\nWith following parameters:\n") ;
 901  //   printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ;
 902  //   printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;
 903  //   printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ;
 904  //   printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;
 905     printf(" ---------------------------------------------------\n") ;
 906   }
 907   else
 908     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
 909
 910 }
 911
 912 //__________________________________________________________________
 913  void AliPHOSDigitizer::PrintDigits(Option_t * option)
 914 {
 915   // Print a table of digits
 916
 917
 918   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 919   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 920   TClonesArray * digits = phosLoader->Digits() ;
 921   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 922
 923   AliInfo(Form("%d", digits->GetEntriesFast())) ;
 924   printf("\nevent %d", gAlice->GetEvNumber()) ;
 925   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;
 926
 927
 928   if(strstr(option,"all")||strstr(option,"EMC")){
 929     //loop over digits
 930     AliPHOSDigit * digit;
 931     printf("\nEMC 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()) &&
 936            (static_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
 937       digit = (AliPHOSDigit * )  digits->At(index) ;
 938       if(digit->GetNprimary() == 0)
 939         continue;
 940 //       printf("%6d  %8d    %6.5e %4d      %2d :",
 941 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
 942       printf("%6d  %.4f    %6.5e %4d      %2d :",
 943               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 944       Int_t iprimary;
 945       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 946         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 947       }
 948       printf("\n") ;
 949     }
 950   }
 951
 952   if(strstr(option,"all")||strstr(option,"CPV")){
 953
 954     //loop over CPV digits
 955     AliPHOSDigit * digit;
 956     printf("\nCPV digits (with primaries):\n")  ;
 957     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 958     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
 959     Int_t index ;
 960     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
 961       digit = (AliPHOSDigit * )  digits->At(index) ;
 962       if(digit->GetId() > maxEmc){
 963         printf("%6d  %8d    %4d      %2d :",
 964                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
 965         Int_t iprimary;
 966         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 967           printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 968         }
 969         printf("\n") ;
 970       }
 971     }
 972   }
 973
 974 }
 975
 976 //__________________________________________________________________
 977 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
 978 {  // Calculates the time signal generated by noise
 979   //PH  Info("TimeOfNoise", "Change me") ;
 980   return gRandom->Rndm() * 1.28E-5;
 981 }
 982
 983 //__________________________________________________________________
 984 void AliPHOSDigitizer::Unload()
 985 {
 986
 987   Int_t i ;
 988   for(i = 1 ; i < fInput ; i++){
 989     TString tempo(fEventNames[i]) ;
 990     tempo += i ;
 991     AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ;
 992     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 993     phosLoader->UnloadSDigits() ;
 994   }
 995
 996   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 997   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 998   phosLoader->UnloadDigits() ;
 999 }
1000
1001 //____________________________________________________________________________
1002 void AliPHOSDigitizer::WriteDigits()
1003 {
1004
1005   // Makes TreeD in the output file.
1006   // Check if branch already exists:
1007   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1008   //      already existing branches.
1009   //   else creates branch with Digits, named "PHOS", title "...",
1010   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
1011   //      and names of files, from which digits are made.
1012
1013   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
1014   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1015
1016   const TClonesArray * digits = phosLoader->Digits() ;
1017   TTree * treeD = phosLoader->TreeD();
1018   if(!treeD){
1019     phosLoader->MakeTree("D");
1020     treeD = phosLoader->TreeD();
1021   }
1022
1023   // -- create Digits branch
1024   Int_t bufferSize = 32000 ;
1025   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
1026   digitsBranch->SetTitle(fEventFolderName);
1027   digitsBranch->Fill() ;
1028
1029   phosLoader->WriteDigits("OVERWRITE");
1030   phosLoader->WriteDigitizer("OVERWRITE");
1031
1032   Unload() ;
1033
1034 }