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 calibration = fcdb->GetADCchannelEmc(module,column,row);
 680     Float_t energy = digit->GetEnergy()/calibration;
 681     digit->SetEnergy(energy); //Now digit measures E in ADC counts
 682     Float_t time = digit->GetTime() ;
 683     time-=fcdb->GetTimeShiftEmc(module,column,row);
 684     digit->SetTime(time) ;
 685   }
 686 }
 687 //____________________________________________________________________________
 688 Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
 689   //Apply time calibration
 690   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 691
 692   //Determine rel.position of the cell absolute ID
 693   Int_t relId[4];
 694   geom->AbsToRelNumbering(absId,relId);
 695   Int_t module=relId[0];
 696   Int_t row   =relId[2];
 697   Int_t column=relId[3];
 698   time += fcdb->GetTimeShiftEmc(module,column,row);
 699   return time ;
 700 }
 701 //____________________________________________________________________________
 702 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
 703 {
 704   // Returns digitized value of the CPV charge in a pad absId
 705
 706   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 707
 708   //Determine rel.position of the cell absId
 709   Int_t relId[4];
 710   geom->AbsToRelNumbering(absId,relId);
 711   Int_t module=relId[0];
 712   Int_t row   =relId[2];
 713   Int_t column=relId[3];
 714
 715   Int_t channel = 0;
 716
 717   if(absId > fEmcCrystals){ //digitize CPV only
 718
 719     //reading calibration data for cell absId.
 720     Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
 721     Float_t adcChanelCpv   = fcdb->GetADCchannelCpv( module,column,row);
 722
 723     channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;
 724     Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
 725     if(channel > nMax ) channel = nMax ;
 726   }
 727   return channel ;
 728 }
 729
 730 //____________________________________________________________________________
 731 void AliPHOSDigitizer::Exec(Option_t *option)
 732 {
 733   // Steering method to process digitization for events
 734   // in the range from fFirstEvent to fLastEvent.
 735   // This range is optionally set by SetEventRange().
 736   // if fLastEvent=-1, then process events until the end.
 737   // by default fLastEvent = fFirstEvent (process only one event)
 738
 739   if (!fInit) { // to prevent overwrite existing file
 740     AliError(Form("Give a version name different from %s",
 741                   fEventFolderName.Data() )) ;
 742     return ;
 743   }
 744
 745   if (strstr(option,"print")) {
 746     Print();
 747     return ;
 748   }
 749
 750  if(strstr(option,"tim"))
 751      gBenchmark->Start("PHOSDigitizer");
 752
 753   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 754   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 755
 756   // Post Digitizer to the white board
 757   phosLoader->PostDigitizer(this) ;
 758
 759   if (fLastEvent == -1)
 760     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 761   else if (fManager)
 762     fLastEvent = fFirstEvent ;
 763
 764   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 765
 766   Int_t ievent ;
 767
 768   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 769     fEventCounter++ ;
 770
 771     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 772
 773     WriteDigits() ;
 774
 775     if(strstr(option,"deb"))
 776       PrintDigits(option);
 777
 778     //increment the total number of Digits per run
 779     fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;
 780  }
 781
 782   phosLoader->CleanDigitizer();
 783
 784   if(strstr(option,"tim")){
 785     gBenchmark->Stop("PHOSDigitizer");
 786     TString message ;
 787     message = "  took %f seconds for Digitizing %f seconds per event\n" ;
 788     AliInfo(Form( message.Data(),
 789          gBenchmark->GetCpuTime("PHOSDigitizer"),
 790          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
 791   }
 792 }
 793 //____________________________________________________________________________
 794 Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
 795   //calculate TOF resolution using beam-test resutls
 796   Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
 797   Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
 798   return TMath::Sqrt(a*a+b*b/e) ;
 799 }
 800
 801 ////____________________________________________________________________________
 802 //Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
 803 //{
 804 //  // Returns the shortest time among all time ticks
 805 //
 806 //  ticks->Sort() ; //Sort in accordance with times of ticks
 807 //  TIter it(ticks) ;
 808 //  AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
 809 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 810 //
 811 //  AliPHOSTick * t ;
 812 //  while((t=(AliPHOSTick*) it.Next())){
 813 //    if(t->GetTime() < time)  //This tick starts before crossing
 814 //      *ctick+=*t ;
 815 //    else
 816 //      return time ;
 817 //
 818 //    time = ctick->CrossingTime(fTimeThreshold) ;
 819 //  }
 820 //  return time ;
 821 //}
 822
 823 //____________________________________________________________________________
 824 Bool_t AliPHOSDigitizer::Init()
 825 {
 826   // Makes all memory allocations
 827   fInit = kTRUE ;
 828
 829   AliPHOSGeometry *geom;
 830   if (!(geom = AliPHOSGeometry::GetInstance()))
 831         geom = AliPHOSGeometry::GetInstance("IHEP","");
 832 //   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 833
 834   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
 835
 836   fFirstEvent = 0 ;
 837   fLastEvent = fFirstEvent ;
 838   if (fManager)
 839     fInput = fManager->GetNinputs() ;
 840   else
 841     fInput           = 1 ;
 842
 843   fInputFileNames  = new TString[fInput] ;
 844   fEventNames      = new TString[fInput] ;
 845   fInputFileNames[0] = GetTitle() ;
 846   fEventNames[0]     = fEventFolderName.Data() ;
 847   Int_t index ;
 848   for (index = 1 ; index < fInput ; index++) {
 849     fInputFileNames[index] = static_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 850     TString tempo = fManager->GetInputFolderName(index) ;
 851     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
 852   }
 853
 854   //to prevent cleaning of this object while GetEvent is called
 855   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 856   if(!rl){
 857     rl = AliRunLoader::Open(GetTitle(), fEventFolderName) ;
 858   }
 859   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 860   phosLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 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 (fManager)
 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   phosLoader->WriteDigitizer("OVERWRITE");
1030
1031   Unload() ;
1032
1033 }