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