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
  19 //_________________________________________________________________________
  20 //*-- Author :  Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  21 //////////////////////////////////////////////////////////////////////////////
  22 // This TTask performs digitization of Summable digits (in the PHOS case it is just
  23 // the sum of contributions from all primary particles into a given cell).
  24 // In addition it performs mixing of summable digits from different events.
  25 // The name of the TTask is also the title of the branch that will contain
  26 // the created SDigits
  27 // The title of the TTAsk is the name of the file that contains the hits from
  28 // which the SDigits are created
  29 //
  30 // For each event two branches are created in TreeD:
  31 //   "PHOS" - list of digits
  32 //   "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
  33 //
  34 // Note, that one can set a title for new digits branch, and repeat digitization with
  35 // another set of parameters.
  36 //
  37 // Use case:
  38 // root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
  39 // root[1] d->ExecuteTask()
  40 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  41 //                       //Digitizes SDigitis in all events found in file galice.root
  42 //
  43 // root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;
  44 //                       // Will read sdigits from galice1.root
  45 // root[3] d1->MixWith("galice2.root")
  46 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  47 //                       // Reads another set of sdigits from galice2.root
  48 // root[3] d1->MixWith("galice3.root")
  49 //                       // Reads another set of sdigits from galice3.root
  50 // root[4] d->ExecuteTask("deb timing")
  51 //                       // Reads SDigits from files galice1.root, galice2.root ....
  52 //                       // mixes them and stores produced Digits in file galice1.root
  53 //                       // deb - prints number of produced digits
  54 //                       // deb all - prints list of produced digits
  55 //                       // timing  - prints time used for digitization
  56 //
  57
  58 // --- ROOT system ---
  59 #include "TTree.h"
  60 #include "TSystem.h"
  61 #include "TBenchmark.h"
  62 #include "TRandom.h"
  63
  64 // --- Standard library ---
  65
  66 // --- AliRoot header files ---
  67
  68 #include "AliRunDigitizer.h"
  69 #include "AliPHOSDigit.h"
  70 #include "AliPHOSGetter.h"
  71 #include "AliPHOSDigitizer.h"
  72 #include "AliPHOSSDigitizer.h"
  73 #include "AliPHOSGeometry.h"
  74 #include "AliPHOSTick.h"
  75
  76 ClassImp(AliPHOSDigitizer)
  77
  78
  79 //____________________________________________________________________________
  80   AliPHOSDigitizer::AliPHOSDigitizer():AliDigitizer("",""),
  81                                        fInput(0),
  82                                        fInputFileNames(0x0),
  83                                        fEventNames(0x0)
  84 {
  85   // ctor
  86   InitParameters() ;
  87   fDefaultInit = kTRUE ;
  88   fManager = 0 ;                     // We work in the standalong mode
  89   fEventFolderName = "" ;
  90 }
  91
  92 //____________________________________________________________________________
  93 AliPHOSDigitizer::AliPHOSDigitizer(const TString alirunFileName, const TString eventFolderName):
  94   AliDigitizer("PHOS"+AliConfig::fgkDigitizerTaskName, alirunFileName),
  95   fInputFileNames(0), fEventNames(0), fEventFolderName(eventFolderName)
  96 {
  97   // ctor
  98
  99   InitParameters() ;
 100   Init() ;
 101   fDefaultInit = kFALSE ;
 102   fManager = 0 ;                     // We work in the standalong mode
 103 }
 104
 105 //____________________________________________________________________________
 106 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d)
 107   : AliDigitizer(d)
 108 {
 109   // copyy ctor
 110
 111   SetName(d.GetName()) ;
 112   SetTitle(d.GetTitle()) ;
 113   fPinNoise = d.fPinNoise ;
 114   fEMCDigitThreshold = d.fEMCDigitThreshold ;
 115   fCPVNoise          = d.fCPVNoise ;
 116   fCPVDigitThreshold = d.fCPVDigitThreshold ;
 117   fTimeResolution    = d.fTimeResolution ;
 118   fTimeThreshold     = d.fTimeThreshold ;
 119   fTimeSignalLength  = d.fTimeSignalLength ;
 120   fADCchanelEmc      = d.fADCchanelEmc ;
 121   fADCpedestalEmc    = d.fADCpedestalEmc ;
 122   fNADCemc           = d.fNADCemc ;
 123   fADCchanelCpv      = d.fADCchanelCpv ;
 124   fADCpedestalCpv    = d.fADCpedestalCpv ;
 125   fNADCcpv           = d.fNADCcpv ;
 126   fEventFolderName   = d.fEventFolderName;
 127 }
 128
 129 //____________________________________________________________________________
 130 AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd):
 131  AliDigitizer(rd,"PHOS"+AliConfig::fgkDigitizerTaskName),
 132  fEventFolderName(0)
 133 {
 134   // ctor
 135   fManager = rd ;
 136   //  SetName(fManager->GetInputFolderName(0)) ;
 137   // take title as name of stream 0
 138   SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
 139   InitParameters() ;
 140   Init() ;
 141   fDefaultInit = kFALSE ;
 142 }
 143
 144 //____________________________________________________________________________
 145   AliPHOSDigitizer::~AliPHOSDigitizer()
 146 {
 147   // dtor
 148   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(),fEventFolderName);
 149   gime->PhosLoader()->CleanDigitizer();
 150   delete [] fInputFileNames ;
 151   delete [] fEventNames ;
 152
 153 }
 154
 155 //____________________________________________________________________________
 156 void AliPHOSDigitizer::Digitize(const Int_t event)
 157 {
 158
 159   // Makes the digitization of the collected summable digits.
 160   //  It first creates the array of all PHOS modules
 161   //  filled with noise (different for EMC, CPV and PPSD) and
 162   //  then adds contributions from SDigits.
 163   // This design avoids scanning over the list of digits to add
 164   // contribution to new SDigits only.
 165
 166   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 167   TClonesArray * digits = gime->Digits() ;
 168   digits->Clear() ;
 169
 170   const AliPHOSGeometry *geom = gime->PHOSGeometry() ;
 171   //Making digits with noise, first EMC
 172   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 173
 174   Int_t nCPV ;
 175   Int_t absID ;
 176
 177   nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * geom->GetNModules() ;
 178
 179   digits->Expand(nCPV) ;
 180
 181   // get first the sdigitizer from the tasks list
 182   if ( !gime->SDigitizer() )
 183     gime->LoadSDigitizer();
 184   AliPHOSSDigitizer * sDigitizer = gime->SDigitizer();
 185
 186   if ( !sDigitizer )
 187     Fatal("Digitize", "SDigitizer with name %s %s not found", fEventFolderName.Data(), GetTitle() ) ;
 188
 189   //take all the inputs to add together and load the SDigits
 190   TObjArray * sdigArray = new TObjArray(fInput) ;
 191   sdigArray->AddAt(gime->SDigits(), 0) ;
 192   Int_t i ;
 193   for(i = 1 ; i < fInput ; i++){
 194     TString tempo(fEventNames[i]) ;
 195     tempo += i ;
 196     AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileNames[i], tempo) ;
 197     gime->Event(event,"S");
 198     sdigArray->AddAt(gime->SDigits(), i) ;
 199   }
 200
 201   //Find the first crystall with signal
 202   Int_t nextSig = 200000 ;
 203   TClonesArray * sdigits ;
 204   for(i = 0 ; i < fInput ; i++){
 205     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 206     if ( !sdigits->GetEntriesFast() )
 207       continue ;
 208     Int_t curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
 209     if(curNext < nextSig)
 210       nextSig = curNext ;
 211   }
 212
 213   TArrayI index(fInput) ;
 214   index.Reset() ;  //Set all indexes to zero
 215
 216   AliPHOSDigit * digit ;
 217   AliPHOSDigit * curSDigit ;
 218
 219   TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
 220
 221   //Put Noise contribution
 222   for(absID = 1 ; absID <= nEMC ; absID++){
 223     Float_t noise = gRandom->Gaus(0., fPinNoise) ;
 224     new((*digits)[absID-1]) AliPHOSDigit( -1, absID, sDigitizer->Digitize(noise), TimeOfNoise() ) ;
 225     //look if we have to add signal?
 226     digit = dynamic_cast<AliPHOSDigit *>(digits->At(absID-1)) ;
 227
 228     if(absID==nextSig){
 229       //Add SDigits from all inputs
 230       ticks->Clear() ;
 231       Int_t contrib = 0 ;
 232       Float_t a = digit->GetAmp() ;
 233       Float_t b = TMath::Abs( a / fTimeSignalLength) ;
 234       //Mark the beginning of the signal
 235       new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);
 236       //Mark the end of the ignal
 237       new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b);
 238
 239       //loop over inputs
 240       for(i = 0 ; i < fInput ; i++){
 241         if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 242           curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 243         else
 244           curSDigit = 0 ;
 245         //May be several digits will contribute from the same input
 246         while(curSDigit && curSDigit->GetId() == absID){
 247           //Shift primary to separate primaries belonging different inputs
 248           Int_t primaryoffset ;
 249           if(fManager)
 250             primaryoffset = fManager->GetMask(i) ;
 251           else
 252             primaryoffset = 10000000*i ;
 253           curSDigit->ShiftPrimary(primaryoffset) ;
 254
 255           a = curSDigit->GetAmp() ;
 256           b = a /fTimeSignalLength ;
 257           new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);
 258           new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
 259
 260           *digit = *digit + *curSDigit ;  //add energies
 261
 262           index[i]++ ;
 263           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 264             curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 265           else
 266             curSDigit = 0 ;
 267         }
 268       }
 269
 270       //calculate and set time
 271       Float_t time = FrontEdgeTime(ticks) ;
 272       digit->SetTime(time) ;
 273
 274       //Find next signal module
 275       nextSig = 200000 ;
 276       for(i = 0 ; i < fInput ; i++){
 277         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 278         Int_t curNext = nextSig ;
 279         if(sdigits->GetEntriesFast() > index[i] ){
 280           curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
 281         }
 282         if(curNext < nextSig) nextSig = curNext ;
 283       }
 284     }
 285   }
 286
 287   ticks->Delete() ;
 288   delete ticks ;
 289
 290   //Now CPV digits (different noise and no timing)
 291   for(absID = nEMC+1; absID <= nCPV; absID++){
 292     Float_t noise = gRandom->Gaus(0., fCPVNoise) ;
 293     new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
 294     //look if we have to add signal?
 295     if(absID==nextSig){
 296       digit = dynamic_cast<AliPHOSDigit *>(digits->At(absID-1)) ;
 297       //Add SDigits from all inputs
 298       for(i = 0 ; i < fInput ; i++){
 299         if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 300           curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 301         else
 302           curSDigit = 0 ;
 303
 304         //May be several digits will contribute from the same input
 305         while(curSDigit && curSDigit->GetId() == absID){
 306           //Shift primary to separate primaries belonging different inputs
 307           Int_t primaryoffset ;
 308           if(fManager)
 309             primaryoffset = fManager->GetMask(i) ;
 310           else
 311             primaryoffset = 10000000*i ;
 312           curSDigit->ShiftPrimary(primaryoffset) ;
 313
 314           //add energies
 315           *digit = *digit + *curSDigit ;
 316           index[i]++ ;
 317           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 318             curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;
 319           else
 320             curSDigit = 0 ;
 321         }
 322       }
 323
 324       //Find next signal module
 325       nextSig = 200000 ;
 326       for(i = 0 ; i < fInput ; i++){
 327         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 328         Int_t curNext = nextSig ;
 329         if(sdigits->GetEntriesFast() > index[i] )
 330           curNext = dynamic_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
 331         if(curNext < nextSig) nextSig = curNext ;
 332       }
 333
 334     }
 335   }
 336
 337   delete sdigArray ; //We should not delete its contents
 338
 339   //remove digits below thresholds
 340   for(i = 0 ; i < nEMC ; i++){
 341     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
 342     if(sDigitizer->Calibrate( digit->GetAmp() ) < fEMCDigitThreshold)
 343       digits->RemoveAt(i) ;
 344     else
 345       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
 346   }
 347
 348
 349   for(i = nEMC; i < nCPV ; i++)
 350     if( sDigitizer->Calibrate( dynamic_cast<AliPHOSDigit*>(digits->At(i))->GetAmp() ) < fCPVDigitThreshold )
 351       digits->RemoveAt(i) ;
 352
 353   digits->Compress() ;
 354
 355   Int_t ndigits = digits->GetEntriesFast() ;
 356   digits->Expand(ndigits) ;
 357
 358   //Set indexes in list of digits and make true digitization of the energy
 359   for (i = 0 ; i < ndigits ; i++) {
 360     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
 361     digit->SetIndexInList(i) ;
 362     Float_t energy = sDigitizer->Calibrate(digit->GetAmp()) ;
 363     digit->SetAmp(DigitizeEnergy(energy,digit->GetId()) ) ;
 364   }
 365 }
 366
 367 //____________________________________________________________________________
 368 Int_t AliPHOSDigitizer::DigitizeEnergy(Float_t energy, Int_t absId)
 369 {
 370   // Returns digitized value of the energy in a cell absId
 371
 372   Int_t chanel ;
 373   if(absId <= fEmcCrystals){ //digitize as EMC
 374     chanel = (Int_t) TMath::Ceil((energy - fADCpedestalEmc)/fADCchanelEmc) ;
 375     if(chanel > fNADCemc ) chanel =  fNADCemc ;
 376   }
 377   else{ //Digitize as CPV
 378     chanel = (Int_t) TMath::Ceil((energy - fADCpedestalCpv)/fADCchanelCpv) ;
 379     if(chanel > fNADCcpv ) chanel =  fNADCcpv ;
 380   }
 381   return chanel ;
 382 }
 383
 384 //____________________________________________________________________________
 385 void AliPHOSDigitizer::Exec(Option_t *option)
 386 {
 387   // Does the job
 388
 389   if (!fInit) { // to prevent overwrite existing file
 390     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
 391     return ;
 392   }
 393
 394   if (strstr(option,"print")) {
 395     Print();
 396     return ;
 397   }
 398
 399   if(strstr(option,"tim"))
 400     gBenchmark->Start("PHOSDigitizer");
 401
 402   if (fManager)
 403     fInput = fManager->GetNinputs() ;
 404
 405   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 406
 407   Int_t nevents = gime->MaxEvent() ;
 408
 409   Int_t ievent ;
 410
 411   for(ievent = 0; ievent < nevents; ievent++){
 412
 413     gime->Event(ievent,"S") ;
 414
 415     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 416
 417     WriteDigits() ;
 418
 419     if(strstr(option,"deb"))
 420       PrintDigits(option);
 421
 422     //increment the total number of Digits per run
 423     fDigitsInRun += gime->Digits()->GetEntriesFast() ;
 424  }
 425
 426   if(strstr(option,"tim")){
 427     gBenchmark->Stop("PHOSDigitizer");
 428     TString message ;
 429     message = "  took %f seconds for Digitizing %f seconds per event\n" ;
 430     Info("Exec", message.Data(),
 431          gBenchmark->GetCpuTime("PHOSDigitizer"),
 432          gBenchmark->GetCpuTime("PHOSDigitizer")/nevents );
 433   }
 434 }
 435
 436 //____________________________________________________________________________
 437 Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
 438 {
 439   // Returns the shortest time among all time ticks
 440
 441   ticks->Sort() ; //Sort in accordance with times of ticks
 442   TIter it(ticks) ;
 443   AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
 444   Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 445
 446   AliPHOSTick * t ;
 447   while((t=(AliPHOSTick*) it.Next())){
 448     if(t->GetTime() < time)  //This tick starts before crossing
 449       *ctick+=*t ;
 450     else
 451       return time ;
 452
 453     time = ctick->CrossingTime(fTimeThreshold) ;
 454   }
 455   return time ;
 456 }
 457
 458 //____________________________________________________________________________
 459 Bool_t AliPHOSDigitizer::Init()
 460 {
 461   // Makes all memory allocations
 462   fInit = kTRUE ;
 463   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 464   if ( gime == 0 ) {
 465     Fatal("Init" ,"Could not obtain the Getter object for file %s and event %s !", GetTitle(), fEventFolderName.Data()) ;
 466     return kFALSE;
 467   }
 468
 469   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 470
 471   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
 472
 473   TString opt("Digits") ;
 474   if(gime->VersionExists(opt) ) {
 475     Error( "Init", "Give a version name different from %s", fEventFolderName.Data() ) ;
 476     fInit = kFALSE ;
 477   }
 478
 479   // Post Digitizer to the white board
 480   gime->PostDigitizer(this) ;
 481
 482   if (fManager)
 483     fInput = fManager->GetNinputs() ;
 484   else
 485     fInput           = 1 ;
 486
 487   fInputFileNames  = new TString[fInput] ;
 488   fEventNames      = new TString[fInput] ;
 489   fInputFileNames[0] = GetTitle() ;
 490   fEventNames[0]     = fEventFolderName.Data() ;
 491   Int_t index ;
 492   for (index = 1 ; index < fInput ; index++) {
 493     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
 494     TString tempo = fManager->GetInputFolderName(index) ;
 495     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager
 496   }
 497
 498   //to prevent cleaning of this object while GetEvent is called
 499   gime->PhosLoader()->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 500
 501   return fInit ;
 502 }
 503
 504 //____________________________________________________________________________
 505 void AliPHOSDigitizer::InitParameters()
 506 {
 507   // Set initial parameters for Digitizer
 508
 509   fPinNoise           = 0.004 ;
 510   fEMCDigitThreshold  = 0.012 ;
 511   fCPVNoise           = 0.01;
 512   fCPVDigitThreshold  = 0.09 ;
 513   fTimeResolution     = 0.5e-9 ;
 514   fTimeSignalLength   = 1.0e-9 ;
 515   fDigitsInRun  = 0 ;
 516   fADCchanelEmc = 0.0015;        // width of one ADC channel in GeV
 517   fADCpedestalEmc = 0.005 ;      //
 518   fNADCemc = (Int_t) TMath::Power(2,16) ;  // number of channels in EMC ADC
 519
 520   fADCchanelCpv = 0.0012 ;          // width of one ADC channel in CPV 'popugais'
 521   fADCpedestalCpv = 0.012 ;         //
 522   fNADCcpv = (Int_t) TMath::Power(2,12);      // number of channels in CPV ADC
 523
 524   fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
 525
 526 }
 527
 528 //__________________________________________________________________
 529 void AliPHOSDigitizer::MixWith(const TString alirunFileName, const TString eventFolderName)
 530 {
 531   // Allows to produce digits by superimposing background and signal event.
 532   // It is assumed, that headers file with SIGNAL events is opened in
 533   // the constructor.
 534   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed
 535   // Thus we avoid writing (changing) huge and expensive
 536   // backgound files: all output will be writen into SIGNAL, i.e.
 537   // opened in constructor file.
 538   //
 539   // One can open as many files to mix with as one needs.
 540   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
 541   // can be mixed.
 542
 543   if( strcmp(fEventFolderName, "") == 0 )
 544     Init() ;
 545
 546   if(fManager){
 547     Warning("MixWith", "Cannot use this method with AliRunDigitizer\n" ) ;
 548     return ;
 549   }
 550   // looking for file which contains AliRun
 551   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
 552     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
 553     return ;
 554   }
 555   // looking for the file which contains SDigits
 556   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 557   TString fileName( gime->GetSDigitsFileName() ) ;
 558     if ( eventFolderName != AliConfig::fgkDefaultEventFolderName) // only if not the default folder name
 559       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
 560     if ( (gSystem->AccessPathName(fileName)) ) {
 561       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
 562       return ;
 563     }
 564     // need to increase the arrays
 565     TString tempo = fInputFileNames[fInput-1] ;
 566     delete [] fInputFileNames ;
 567     fInputFileNames = new TString[fInput+1] ;
 568     fInputFileNames[fInput-1] = tempo ;
 569
 570     tempo = fEventNames[fInput-1] ;
 571     delete [] fEventNames ;
 572     fEventNames = new TString[fInput+1] ;
 573     fEventNames[fInput-1] = tempo ;
 574
 575     fInputFileNames[fInput] = alirunFileName ;
 576     fEventNames[fInput]     = eventFolderName ;
 577     fInput++ ;
 578 }
 579
 580 //__________________________________________________________________
 581 void AliPHOSDigitizer::Print()const
 582 {
 583   // Print Digitizer's parameters
 584   Info("Print", "\n------------------- %s -------------", GetName() ) ;
 585   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 586     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 587
 588     Int_t nStreams ;
 589     if (fManager)
 590       nStreams =  GetNInputStreams() ;
 591     else
 592       nStreams = fInput ;
 593
 594     Int_t index = 0 ;
 595     for (index = 0 ; index < nStreams ; index++) {
 596       TString tempo(fEventNames[index]) ;
 597       tempo += index ;
 598       AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileNames[index], tempo) ;
 599       TString fileName( gime->GetSDigitsFileName() ) ;
 600       if ( fEventNames[index] != AliConfig::fgkDefaultEventFolderName) // only if not the default folder name
 601         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
 602       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
 603     }
 604     AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 605     printf("\nWriting digits to %s", gime->GetDigitsFileName().Data()) ;
 606
 607     printf("\nWith following parameters:\n") ;
 608     printf("  Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise ) ;
 609     printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f\n", fEMCDigitThreshold ) ;
 610     printf("  Noise in CPV (fCPVNoise) = %f\n", fCPVNoise ) ;
 611     printf("  Threshold in CPV (fCPVDigitThreshold) = %f\n",fCPVDigitThreshold ) ;
 612     printf(" ---------------------------------------------------\n") ;
 613   }
 614   else
 615     Info("Print", "AliPHOSDigitizer not initialized" ) ;
 616
 617 }
 618
 619 //__________________________________________________________________
 620  void AliPHOSDigitizer::PrintDigits(Option_t * option)
 621 {
 622   // Print a table of digits
 623
 624   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 625   TClonesArray * digits = gime->Digits() ;
 626
 627   Info("PrintDigits", "%d", digits->GetEntriesFast()) ;
 628   printf("\nevent %d", gAlice->GetEvNumber()) ;
 629   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;
 630
 631
 632   if(strstr(option,"all")||strstr(option,"EMC")){
 633     //loop over digits
 634     AliPHOSDigit * digit;
 635     printf("\nEMC digits (with primaries):\n")  ;
 636     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 637     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
 638     Int_t index ;
 639     for (index = 0 ; (index < digits->GetEntriesFast()) &&
 640            (dynamic_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
 641       digit = (AliPHOSDigit * )  digits->At(index) ;
 642       if(digit->GetNprimary() == 0)
 643         continue;
 644       printf("%6d  %8d    %6.5e %4d      %2d :",
 645               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 646       Int_t iprimary;
 647       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 648         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 649       }
 650       printf("\n") ;
 651     }
 652   }
 653
 654   if(strstr(option,"all")||strstr(option,"CPV")){
 655
 656     //loop over CPV digits
 657     AliPHOSDigit * digit;
 658     printf("\nCPV digits (with primaries):\n")  ;
 659     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 660     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
 661     Int_t index ;
 662     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
 663       digit = (AliPHOSDigit * )  digits->At(index) ;
 664       if(digit->GetId() > maxEmc){
 665         printf("%6d  %8d    %4d      %2d :",
 666                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
 667         Int_t iprimary;
 668         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 669           printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 670         }
 671         printf("\n") ;
 672       }
 673     }
 674   }
 675
 676 }
 677
 678 //__________________________________________________________________
 679 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
 680 {  // Calculates the time signal generated by noise
 681   //to be rewritten, now returns just big number
 682   return 1. ;
 683
 684 }
 685
 686 //__________________________________________________________________
 687 void AliPHOSDigitizer::Unload()
 688 {
 689
 690   Int_t i ;
 691   for(i = 1 ; i < fInput ; i++){
 692     TString tempo(fEventNames[i]) ;
 693     tempo += i ;
 694     AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileNames[i], tempo) ;
 695     gime->PhosLoader()->UnloadSDigits() ;
 696   }
 697
 698   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 699   gime->PhosLoader()->UnloadDigits() ;
 700 }
 701
 702 //____________________________________________________________________________
 703 void AliPHOSDigitizer::WriteDigits()
 704 {
 705
 706   // Makes TreeD in the output file.
 707   // Check if branch already exists:
 708   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
 709   //      already existing branches.
 710   //   else creates branch with Digits, named "PHOS", title "...",
 711   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
 712   //      and names of files, from which digits are made.
 713
 714   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ;
 715   const TClonesArray * digits = gime->Digits() ;
 716   TTree * treeD = gime->TreeD();
 717
 718   // -- create Digits branch
 719   Int_t bufferSize = 32000 ;
 720   TBranch * digitsBranch = treeD->Branch("PHOS",&digits,bufferSize);
 721   digitsBranch->SetTitle(fEventFolderName);
 722   digitsBranch->Fill() ;
 723
 724   gime->WriteDigits("OVERWRITE");
 725   gime->WriteDigitizer("OVERWRITE");
 726
 727   Unload() ;
 728
 729 }