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