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