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