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