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