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