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