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