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: AliPHOSDigitizer.cxx,v $
  21  * Revision 1.104  2007/12/18 09:08:18  hristov
  22  * Splitting of the QA maker into simulation and reconstruction dependent parts (Yves)
  23  *
  24  * Revision 1.103  2007/11/07 11:25:06  schutz
  25  * Comment out the QA checking before starting digitization
  26  *
  27  * Revision 1.102  2007/10/19 18:04:29  schutz
  28  * The standalone QA data maker is called from AliSimulation and AliReconstruction outside the event loop; i.e. re-reading the data. The QA data making in the event loop has been commented out.
  29  *
  30  * Revision 1.101  2007/10/14 21:08:10  schutz
  31  * Introduced the checking of QA results from previous step before entering the event loop
  32  *
  33  * Revision 1.100  2007/10/10 09:05:10  schutz
  34  * Changing name QualAss to QA
  35  *
  36  * Revision 1.99  2007/09/30 17:08:20  schutz
  37  * Introducing the notion of QA data acquisition cycle (needed by online)
  38  *
  39  * Revision 1.98  2007/09/26 14:22:17  cvetan
  40  * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
  41  *
  42  * Revision 1.97  2007/08/07 14:12:03  kharlov
  43  * Quality assurance added (Yves Schutz)
  44  *
  45  * Revision 1.96  2007/04/28 10:43:36  policheh
  46  * Dead channels simulation: digit energy sets to 0.
  47  *
  48  * Revision 1.95  2007/04/10 07:20:52  kharlov
  49  * Decalibration should use the same CDB as calibration in AliPHOSClusterizerv1
  50  *
  51  * Revision 1.94  2007/02/01 10:34:47  hristov
  52  * Removing warnings on Solaris x86
  53  *
  54  * Revision 1.93  2006/10/17 13:17:01  kharlov
  55  * Replace AliInfo by AliDebug
  56  *
  57  * Revision 1.92  2006/08/28 10:01:56  kharlov
  58  * Effective C++ warnings fixed (Timur Pocheptsov)
  59  *
  60  * Revision 1.91  2006/04/29 20:25:30  hristov
  61  * Decalibration is implemented (Yu.Kharlov)
  62  *
  63  * Revision 1.90  2006/04/22 10:30:17  hristov
  64  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
  65  *
  66  * Revision 1.89  2006/04/11 15:22:59  hristov
  67  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
  68  *
  69  * Revision 1.88  2006/03/13 14:05:43  kharlov
  70  * Calibration objects for EMC and CPV
  71  *
  72  * Revision 1.87  2005/08/24 15:33:49  kharlov
  73  * Calibration data for raw digits
  74  *
  75  * Revision 1.86  2005/07/12 20:07:35  hristov
  76  * Changes needed to run simulation and reconstrruction in the same AliRoot session
  77  *
  78  * Revision 1.85  2005/05/28 14:19:04  schutz
  79  * Compilation warnings fixed by T.P.
  80  *
  81  */
  82
  83 //_________________________________________________________________________
  84 //*-- Author :  Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  85 //////////////////////////////////////////////////////////////////////////////
  86 // This class performs digitization of Summable digits (in the PHOS case it is just
  87 // the sum of contributions from all primary particles into a given cell).
  88 // In addition it performs mixing of summable digits from different events.
  89 // The name of the class is also the title of the branch that will contain
  90 // the created SDigits
  91 // The title of the class is the name of the file that contains the hits from
  92 // which the SDigits are created
  93 //
  94 // For each event two branches are created in TreeD:
  95 //   "PHOS" - list of digits
  96 //   "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
  97 //
  98 // Note, that one can set a title for new digits branch, and repeat digitization with
  99 // another set of parameters.
 100 //
 101 // Use case:
 102 // root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
 103 // root[1] d->Digitize()
 104 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
 105 //                       //Digitizes SDigitis in all events found in file galice.root
 106 //
 107 // root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;
 108 //                       // Will read sdigits from galice1.root
 109 // root[3] d1->MixWith("galice2.root")
 110 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
 111 //                       // Reads another set of sdigits from galice2.root
 112 // root[3] d1->MixWith("galice3.root")
 113 //                       // Reads another set of sdigits from galice3.root
 114 // root[4] d->Digitize("deb timing")
 115 //                       // Reads SDigits from files galice1.root, galice2.root ....
 116 //                       // mixes them and stores produced Digits in file galice1.root
 117 //                       // deb - prints number of produced digits
 118 //                       // deb all - prints list of produced digits
 119 //                       // timing  - prints time used for digitization
 120 //
 121
 122 // --- ROOT system ---
 123 #include "TTree.h"
 124 #include "TSystem.h"
 125 #include "TBenchmark.h"
 126 #include "TRandom.h"
 127 #include "TMath.h"
 128
 129 // --- Standard library ---
 130
 131 // --- AliRoot header files ---
 132 #include <TGeoManager.h>
 133 #include "AliLog.h"
 134 #include "AliDigitizationInput.h"
 135 #include "AliPHOSDigit.h"
 136 #include "AliPHOSDigitizer.h"
 137 #include "AliPHOSGeometry.h"
 138 #include "AliPHOSTick.h"
 139 #include "AliPHOSSimParam.h"
 140 #include "AliPHOSCalibData.h"
 141 #include "AliRunLoader.h"
 142 #include "AliPHOSLoader.h"
 143 #include "AliPHOSPulseGenerator.h"
 144
 145 ClassImp(AliPHOSDigitizer)
 146
 147
 148 //____________________________________________________________________________
 149 AliPHOSDigitizer::AliPHOSDigitizer() :
 150   AliDigitizer("",""),
 151   fDefaultInit(kTRUE),
 152   fDigitsInRun(0),
 153   fInit(kFALSE),
 154   fInput(0),
 155   fInputFileNames(0x0),
 156   fEventNames(0x0),
 157   fEmcCrystals(0),
 158   fEventFolderName(""),
 159   fFirstEvent(0),
 160   fLastEvent(0),
 161   fcdb(0x0),
 162   fEventCounter(0),
 163   fPulse(0),
 164   fADCValuesLG(0),
 165   fADCValuesHG(0)
 166 {
 167   // ctor
 168   InitParameters() ;
 169   fDigInput = 0 ;                     // We work in the standalong mode
 170 }
 171
 172 //____________________________________________________________________________
 173 AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName,
 174                                    TString eventFolderName):
 175   AliDigitizer("PHOSDigitizer", alirunFileName),
 176   fDefaultInit(kFALSE),
 177   fDigitsInRun(0),
 178   fInit(kFALSE),
 179   fInput(0),
 180   fInputFileNames(0x0),
 181   fEventNames(0x0),
 182   fEmcCrystals(0),
 183   fEventFolderName(eventFolderName),
 184   fFirstEvent(0),
 185   fLastEvent(0),
 186   fcdb(0x0),
 187   fEventCounter(0),
 188   fPulse(0),
 189   fADCValuesLG(0),
 190   fADCValuesHG(0)
 191 {
 192   // ctor
 193   InitParameters() ;
 194   Init() ;
 195   fDefaultInit = kFALSE ;
 196   fDigInput = 0 ;                     // We work in the standalone mode
 197   fcdb = new AliPHOSCalibData(-1);
 198 }
 199
 200 //____________________________________________________________________________
 201 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) :
 202   AliDigitizer(d),
 203   fDefaultInit(d.fDefaultInit),
 204   fDigitsInRun(d.fDigitsInRun),
 205   fInit(d.fInit),
 206   fInput(d.fInput),
 207   fInputFileNames(0x0),
 208   fEventNames(0x0),
 209   fEmcCrystals(d.fEmcCrystals),
 210   fEventFolderName(d.fEventFolderName),
 211   fFirstEvent(d.fFirstEvent),
 212   fLastEvent(d.fLastEvent),
 213   fcdb (new AliPHOSCalibData(-1)),
 214   fEventCounter(0),
 215   fPulse(new AliPHOSPulseGenerator),
 216   fADCValuesLG(0),
 217   fADCValuesHG(0)
 218 {
 219   // copyy ctor
 220   SetName(d.GetName()) ;
 221   SetTitle(d.GetTitle()) ;
 222   for (Int_t iInput=0; iInput<fInput; iInput++) {
 223     fInputFileNames[iInput] = d.fInputFileNames[iInput];
 224     fEventNames[iInput]     = d.fEventNames[iInput];
 225   }
 226 }
 227
 228 //____________________________________________________________________________
 229 AliPHOSDigitizer::AliPHOSDigitizer(AliDigitizationInput * rd) :
 230   AliDigitizer(rd,"PHOSDigitizer"),
 231   fDefaultInit(kFALSE),
 232   fDigitsInRun(0),
 233   fInit(kFALSE),
 234   fInput(0),
 235   fInputFileNames(0x0),
 236   fEventNames(0x0),
 237   fEmcCrystals(0),
 238   fEventFolderName(fDigInput->GetInputFolderName(0)),
 239   fFirstEvent(0),
 240   fLastEvent(0),
 241   fcdb (0x0),
 242   fEventCounter(0),
 243   fPulse(0),
 244   fADCValuesLG(0),
 245   fADCValuesHG(0)
 246
 247 {
 248   // ctor Init() is called by RunDigitizer
 249   fDigInput = rd ;
 250   SetTitle(static_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
 251   InitParameters() ;
 252   fDefaultInit = kFALSE ;
 253   fcdb = new AliPHOSCalibData(-1);
 254 }
 255
 256 //____________________________________________________________________________
 257   AliPHOSDigitizer::~AliPHOSDigitizer()
 258 {
 259   // dtor
 260   delete [] fInputFileNames ;
 261   delete [] fEventNames ;
 262
 263   delete fPulse;
 264   delete [] fADCValuesLG;
 265   delete [] fADCValuesHG;
 266
 267   if(fcdb){ delete fcdb ; fcdb=0;}
 268
 269 }
 270
 271 //____________________________________________________________________________
 272 void AliPHOSDigitizer::Digitize(Int_t event)
 273 {
 274
 275   // Makes the digitization of the collected summable digits.
 276   //  It first creates the array of all PHOS modules
 277   //  filled with noise (different for EMC, and CPV) and
 278   //  then adds contributions from SDigits.
 279   // This design avoids scanning over the list of digits to add
 280   // contribution to new SDigits only.
 281
 282   Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
 283
 284   //First stream
 285   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 286   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 287
 288   Int_t readEvent = event ;
 289   if (fDigInput)
 290     readEvent = static_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
 291   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 292                   readEvent, GetTitle(), fEventFolderName.Data())) ;
 293   rl->GetEvent(readEvent) ;
 294   phosLoader->CleanSDigits() ;
 295   phosLoader->LoadSDigits("READ") ;
 296
 297   //Prepare Output
 298   TClonesArray * digits = phosLoader->Digits() ;
 299   if( !digits ) {
 300     phosLoader->MakeDigitsArray() ;
 301     digits = phosLoader->Digits() ;
 302   }
 303
 304   digits->Clear() ;
 305
 306   //
 307   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 308   //Making digits with noise, first EMC
 309   //Check which PHOS modules are present
 310   Bool_t isPresent[5] ;
 311   TString volpath ;
 312   Int_t nmod=0 ;
 313   for(Int_t i=0; i<5; i++){
 314     volpath = "/ALIC_1/PHOS_";
 315     volpath += i+1;
 316     if (gGeoManager->CheckPath(volpath.Data())) {
 317       isPresent[i]=1 ;
 318       nmod++ ;
 319     }
 320     else{
 321       isPresent[i]=0 ;
 322     }
 323   }
 324
 325   Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
 326
 327   Int_t nCPV ;
 328   Int_t absID ;
 329
 330   //check if CPV exists
 331   Bool_t isCPVpresent=0 ;
 332   for(Int_t i=1; i<=5 && !isCPVpresent; i++){
 333     volpath = "/ALIC_1/PHOS_";
 334     volpath += i;
 335     volpath += "/PCPV_1";
 336     if (gGeoManager->CheckPath(volpath.Data()))
 337       isCPVpresent=1 ;
 338   }
 339
 340   if(isCPVpresent){
 341     nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
 342   }
 343   else{
 344      nCPV = nEMC ;
 345   }
 346
 347   digits->Expand(nCPV) ;
 348
 349   //take all the inputs to add together and load the SDigits
 350   TObjArray * sdigArray = new TObjArray(fInput) ;
 351   sdigArray->AddAt(phosLoader->SDigits(), 0) ;
 352
 353   for(Int_t i = 1 ; i < fInput ; i++){
 354     TString tempo(fEventNames[i]) ;
 355     tempo += i ;
 356     AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
 357     if(!rl2){
 358       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
 359       if (!rl2) {
 360         Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
 361         return ;
 362       }
 363       rl2->LoadHeader();
 364     }
 365     AliPHOSLoader * phosLoader2 = static_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
 366
 367     if(fDigInput){
 368       readEvent = static_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
 369     }
 370     TClonesArray * digs ;
 371     if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
 372       AliInfo(Form("Adding event %d from input stream %d %s %s",
 373                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
 374        rl2->GetEvent(readEvent) ;
 375        phosLoader2->CleanDigits() ;
 376        phosLoader2->LoadDigits("READ") ;
 377        digs = phosLoader2->Digits() ;
 378        toMakeNoise=0 ; //Do not add noise, it is already in stream
 379     }
 380     else{
 381       AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
 382                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
 383        rl2->GetEvent(readEvent) ;
 384        phosLoader2->CleanSDigits() ;
 385        phosLoader2->LoadSDigits("READ") ;
 386        digs = phosLoader2->SDigits() ;
 387     }
 388     sdigArray->AddAt(digs, i) ;
 389   }
 390
 391   //Find the first crystal with signal
 392   Int_t nextSig = 200000 ;
 393   TClonesArray * sdigits ;
 394   for(Int_t i = 0 ; i < fInput ; i++){
 395     sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 396     if ( !sdigits->GetEntriesFast() )
 397       continue ;
 398     Int_t curNext = static_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
 399     if(curNext < nextSig)
 400       nextSig = curNext ;
 401   }
 402
 403   TArrayI index(fInput) ;
 404   index.Reset() ;  //Set all indexes to zero
 405
 406   AliPHOSDigit * digit ;
 407   AliPHOSDigit * curSDigit ;
 408
 409 //  TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
 410
 411   //Put Noise contribution
 412   Double_t apdNoise = 0. ;
 413   if(toMakeNoise)
 414      apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ;
 415
 416   Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
 417   Int_t idigit= 0;
 418   for(Int_t imod=0; imod<5; imod++){
 419     if(!isPresent[imod])
 420       continue ;
 421     Int_t firstAbsId=imod*emcpermod+1 ;
 422     Int_t lastAbsId =(imod+1)*emcpermod ;
 423     for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
 424       Float_t noise = gRandom->Gaus(0.,apdNoise) ;
 425       new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
 426       //look if we have to add signal?
 427       digit = static_cast<AliPHOSDigit *>(digits->At(idigit)) ;
 428       idigit++ ;
 429
 430       if(absID==nextSig){
 431       //Add SDigits from all inputs
 432 //      ticks->Clear() ;
 433 //      Int_t contrib = 0 ;
 434
 435 //New Timing model is necessary
 436 //      Float_t a = digit->GetEnergy() ;
 437 //      Float_t b = TMath::Abs( a / fTimeSignalLength) ;
 438 //      //Mark the beginning of the signal
 439 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);
 440 //      //Mark the end of the signal
 441 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b);
 442
 443 // Calculate time as time of the largest digit
 444         Float_t time = digit->GetTime() ;
 445         Float_t eTime= digit->GetEnergy() ;
 446
 447         //loop over inputs
 448         for(Int_t i = 0 ; i < fInput ; i++){
 449           if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
 450             curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 451             if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
 452               curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
 453               curSDigit->SetTime(CalibrateT(curSDigit->GetTime(),curSDigit->GetId())) ;
 454             }
 455           }
 456           else
 457             curSDigit = 0 ;
 458           //May be several digits will contribute from the same input
 459           while(curSDigit && curSDigit->GetId() == absID){
 460             //Shift primary to separate primaries belonging different inputs
 461             Int_t primaryoffset ;
 462             if(fDigInput)
 463               primaryoffset = fDigInput->GetMask(i) ;
 464             else
 465               primaryoffset = 10000000*i ;
 466             curSDigit->ShiftPrimary(primaryoffset) ;
 467
 468 //New Timing model is necessary
 469 //        a = curSDigit->GetEnergy() ;
 470 //        b = a /fTimeSignalLength ;
 471 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);
 472 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
 473           if(curSDigit->GetEnergy()>eTime){
 474             eTime=curSDigit->GetEnergy() ;
 475             time=curSDigit->GetTime() ;
 476           }
 477             *digit += *curSDigit ;  //add energies
 478
 479             index[i]++ ;
 480             if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 481               curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 482             else
 483               curSDigit = 0 ;
 484           }
 485         }
 486
 487         //calculate and set time
 488 //New Timing model is necessary
 489 //      Float_t time = FrontEdgeTime(ticks) ;
 490         digit->SetTime(time) ;
 491
 492         //Find next signal module
 493         nextSig = 200000 ;
 494         for(Int_t i = 0 ; i < fInput ; i++){
 495           sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 496           Int_t curNext = nextSig ;
 497           if(sdigits->GetEntriesFast() > index[i] ){
 498             curNext = static_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
 499           }
 500           if(curNext < nextSig) nextSig = curNext ;
 501         }
 502       }
 503     }
 504   }
 505
 506
 507   //Apply non-linearity
 508   if(AliPHOSSimParam::GetInstance()->IsCellNonlinearityOn()){ //Apply non-lineairyt on cell level
 509     const Double_t aNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyA() ;
 510     const Double_t bNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyB() ;
 511     const Double_t cNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyC() ;
 512     for(Int_t i = 0 ; i < nEMC ; i++){
 513       digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 514       Double_t e= digit->GetEnergy() ;
 515       // version(1)      digit->SetEnergy(e*(1+a*TMath::Exp(-e/b))) ;
 516       digit->SetEnergy(e*cNL*(1.+aNL*TMath::Exp(-e*e/2./bNL/bNL))) ; //Better agreement with data...
 517     }
 518   }
 519
 520
 521   //distretize energy if necessary
 522   if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
 523     Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
 524     for(Int_t i = 0 ; i < nEMC ; i++){
 525       digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 526       digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
 527     }
 528   }
 529
 530   //Apply decalibration if necessary
 531   for(Int_t i = 0 ; i < nEMC ; i++){
 532     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 533     Decalibrate(digit) ;
 534   }
 535
 536 //  ticks->Delete() ;
 537 //  delete ticks ;
 538
 539   //Now CPV digits (different noise and no timing)
 540   Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
 541   Int_t nEMCtotal=emcpermod*5 ;
 542   Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
 543   if(isCPVpresent){  //CPV is present in current geometry
 544     for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
 545       if(!isPresent[imod])
 546         continue ;
 547       Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
 548       Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
 549       for(absID = firstAbsId; absID <= lastAbsId; absID++){
 550         Float_t noise = gRandom->Gaus(0., cpvNoise) ;
 551         new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
 552         idigit++ ;
 553         //look if we have to add signal?
 554         if(absID==nextSig){
 555           digit = static_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
 556           //Add SDigits from all inputs
 557           for(Int_t i = 0 ; i < fInput ; i++){
 558              if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 559                curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
 560              else
 561                curSDigit = 0 ;
 562
 563              //May be several digits will contribute from the same input
 564              while(curSDigit && curSDigit->GetId() == absID){
 565                //Shift primary to separate primaries belonging different inputs
 566                Int_t primaryoffset ;
 567                if(fDigInput)
 568                  primaryoffset = fDigInput->GetMask(i) ;
 569                else
 570                  primaryoffset = 10000000*i ;
 571                curSDigit->ShiftPrimary(primaryoffset) ;
 572
 573                //add energies
 574                *digit += *curSDigit ;
 575                index[i]++ ;
 576                if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
 577                  curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;
 578                else
 579                  curSDigit = 0 ;
 580              }
 581           }
 582
 583           //Find next signal module
 584           nextSig = 200000 ;
 585           for(Int_t i = 0 ; i < fInput ; i++){
 586             sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
 587             Int_t curNext = nextSig ;
 588             if(sdigits->GetEntriesFast() > index[i] )
 589               curNext = static_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
 590             if(curNext < nextSig) nextSig = curNext ;
 591           }
 592
 593         }
 594       }
 595     }
 596   }
 597
 598   delete sdigArray ; //We should not delete its contents
 599
 600   Int_t relId[4];
 601
 602   //set amplitudes in bad channels to zero
 603
 604   for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
 605     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 606     geom->AbsToRelNumbering(digit->GetId(),relId);
 607     if(relId[1] == 0) // Emc
 608       if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.);
 609   }
 610
 611   //remove digits below thresholds
 612   Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
 613   for(Int_t i = 0 ; i < nEMC ; i++){
 614     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 615
 616     if(digit->GetEnergy() < emcThreshold){
 617       digits->RemoveAt(i) ;
 618       continue ;
 619     }
 620
 621     geom->AbsToRelNumbering(digit->GetId(),relId);
 622
 623 //    digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
 624
 625     Float_t tres = TimeResolution(digit->GetEnergy()) ;
 626     digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
 627
 628     fPulse->Reset();
 629     fPulse->SetAmplitude(digit->GetEnergy()/
 630                          fcdb->GetADCchannelEmc(relId[0],relId[3],relId[2]));
 631     fPulse->SetTZero(digit->GetTimeR());
 632     fPulse->MakeSamples();
 633     fPulse->GetSamples(fADCValuesHG, fADCValuesLG) ;
 634     Int_t nSamples = fPulse->GetRawFormatTimeBins();
 635     digit->SetALTROSamplesHG(nSamples,fADCValuesHG);
 636     digit->SetALTROSamplesLG(nSamples,fADCValuesLG);
 637   }
 638
 639   Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
 640   for(Int_t i = nEMC; i < nCPV ; i++){
 641     if( static_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
 642       digits->RemoveAt(i) ;
 643   }
 644
 645   digits->Compress() ;
 646   Int_t ndigits = digits->GetEntriesFast() ;
 647
 648   //Set indexes in list of digits and make true digitization of the energy
 649   for (Int_t i = 0 ; i < ndigits ; i++) {
 650     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 651     digit->SetIndexInList(i) ;
 652     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
 653       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
 654     }
 655   }
 656
 657 }
 658 //____________________________________________________________________________
 659 Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
 660   //Apply calibration
 661   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 662
 663   //Determine rel.position of the cell absolute ID
 664   Int_t relId[4];
 665   geom->AbsToRelNumbering(absId,relId);
 666   Int_t module=relId[0];
 667   Int_t row   =relId[2];
 668   Int_t column=relId[3];
 669   if(relId[1]==0){ //This Is EMC
 670     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
 671     return amp*calibration ;
 672   }
 673   return 0 ;
 674 }
 675 //____________________________________________________________________________
 676 void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
 677 {
 678   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
 679
 680   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 681
 682   //Determine rel.position of the cell absolute ID
 683   Int_t relId[4];
 684   geom->AbsToRelNumbering(digit->GetId(),relId);
 685   Int_t module=relId[0];
 686   Int_t row   =relId[2];
 687   Int_t column=relId[3];
 688   if(relId[1]==0){ //This Is EMC
 689     Float_t decalib     = fcdb->GetADCchannelEmcDecalib(module,column,row); // O(1)
 690     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row)*decalib;
 691     Float_t energy = digit->GetEnergy()/calibration;
 692     digit->SetEnergy(energy); //Now digit measures E in ADC counts
 693     Float_t time = digit->GetTime() ;
 694     time-=fcdb->GetTimeShiftEmc(module,column,row);
 695     digit->SetTime(time) ;
 696   }
 697 }
 698 //____________________________________________________________________________
 699 Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
 700   //Apply time calibration
 701   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 702
 703   //Determine rel.position of the cell absolute ID
 704   Int_t relId[4];
 705   geom->AbsToRelNumbering(absId,relId);
 706   Int_t module=relId[0];
 707   Int_t row   =relId[2];
 708   Int_t column=relId[3];
 709   time += fcdb->GetTimeShiftEmc(module,column,row);
 710   return time ;
 711 }
 712 //____________________________________________________________________________
 713 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
 714 {
 715   // Returns digitized value of the CPV charge in a pad absId
 716
 717   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 718
 719   //Determine rel.position of the cell absId
 720   Int_t relId[4];
 721   geom->AbsToRelNumbering(absId,relId);
 722   Int_t module=relId[0];
 723   Int_t row   =relId[2];
 724   Int_t column=relId[3];
 725
 726   Int_t channel = 0;
 727
 728   if(absId > fEmcCrystals){ //digitize CPV only
 729
 730     //reading calibration data for cell absId.
 731     Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
 732     Float_t adcChanelCpv   = fcdb->GetADCchannelCpv( module,column,row);
 733
 734     channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;
 735     Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
 736     if(channel > nMax ) channel = nMax ;
 737   }
 738   return channel ;
 739 }
 740
 741 //____________________________________________________________________________
 742 void AliPHOSDigitizer::Digitize(Option_t *option)
 743 {
 744   // Steering method to process digitization for events
 745   // in the range from fFirstEvent to fLastEvent.
 746   // This range is optionally set by SetEventRange().
 747   // if fLastEvent=-1, then process events until the end.
 748   // by default fLastEvent = fFirstEvent (process only one event)
 749
 750   if (!fInit) { // to prevent overwrite existing file
 751     AliError(Form("Give a version name different from %s",
 752                   fEventFolderName.Data() )) ;
 753     return ;
 754   }
 755
 756   if (strstr(option,"print")) {
 757     Print();
 758     return ;
 759   }
 760
 761  if(strstr(option,"tim"))
 762      gBenchmark->Start("PHOSDigitizer");
 763
 764   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 765   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 766
 767   if (fLastEvent == -1)
 768     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 769   else if (fDigInput)
 770     fLastEvent = fFirstEvent ;
 771
 772   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 773
 774   Int_t ievent ;
 775
 776   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 777     fEventCounter++ ;
 778
 779     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 780
 781     WriteDigits() ;
 782
 783     if(strstr(option,"deb"))
 784       PrintDigits(option);
 785
 786     //increment the total number of Digits per run
 787     fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;
 788  }
 789
 790   if(strstr(option,"tim")){
 791     gBenchmark->Stop("PHOSDigitizer");
 792     TString message ;
 793     message = "  took %f seconds for Digitizing %f seconds per event\n" ;
 794     AliInfo(Form( message.Data(),
 795          gBenchmark->GetCpuTime("PHOSDigitizer"),
 796          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
 797   }
 798 }
 799 //____________________________________________________________________________
 800 Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
 801   //calculate TOF resolution using beam-test resutls
 802   Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
 803   Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
 804   return TMath::Sqrt(a*a+b*b/e) ;
 805 }
 806
 807 ////____________________________________________________________________________
 808 //Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
 809 //{
 810 //  // Returns the shortest time among all time ticks
 811 //
 812 //  ticks->Sort() ; //Sort in accordance with times of ticks
 813 //  TIter it(ticks) ;
 814 //  AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
 815 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;
 816 //
 817 //  AliPHOSTick * t ;
 818 //  while((t=(AliPHOSTick*) it.Next())){
 819 //    if(t->GetTime() < time)  //This tick starts before crossing
 820 //      *ctick+=*t ;
 821 //    else
 822 //      return time ;
 823 //
 824 //    time = ctick->CrossingTime(fTimeThreshold) ;
 825 //  }
 826 //  return time ;
 827 //}
 828
 829 //____________________________________________________________________________
 830 Bool_t AliPHOSDigitizer::Init()
 831 {
 832   // Makes all memory allocations
 833   fInit = kTRUE ;
 834
 835   AliPHOSGeometry *geom;
 836   if (!(geom = AliPHOSGeometry::GetInstance()))
 837         geom = AliPHOSGeometry::GetInstance("IHEP","");
 838 //   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 839
 840   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
 841
 842   fFirstEvent = 0 ;
 843   fLastEvent = fFirstEvent ;
 844   if (fDigInput)
 845     fInput = fDigInput->GetNinputs() ;
 846   else
 847     fInput           = 1 ;
 848
 849   fInputFileNames  = new TString[fInput] ;
 850   fEventNames      = new TString[fInput] ;
 851   fInputFileNames[0] = GetTitle() ;
 852   fEventNames[0]     = fEventFolderName.Data() ;
 853   Int_t index ;
 854   for (index = 1 ; index < fInput ; index++) {
 855     fInputFileNames[index] = static_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
 856     TString tempo = fDigInput->GetInputFolderName(index) ;
 857     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fDigInput
 858   }
 859
 860   //to prevent cleaning of this object while GetEvent is called
 861   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 862   if(!rl){
 863     AliRunLoader::Open(GetTitle(), fEventFolderName) ;
 864   }
 865   return fInit ;
 866 }
 867
 868 //____________________________________________________________________________
 869 void AliPHOSDigitizer::InitParameters()
 870 {
 871   // Set initial parameters Digitizer
 872
 873   fDigitsInRun  = 0 ;
 874   SetEventRange(0,-1) ;
 875   fPulse = new AliPHOSPulseGenerator();
 876   fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
 877   fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
 878
 879 }
 880
 881 //__________________________________________________________________
 882 void AliPHOSDigitizer::Print(const Option_t *)const
 883 {
 884   // Print Digitizer's parameters
 885   AliInfo(Form("\n------------------- %s -------------", GetName() )) ;
 886   if( strcmp(fEventFolderName.Data(), "") != 0 ){
 887     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
 888
 889     Int_t nStreams ;
 890     if (fDigInput)
 891       nStreams =  GetNInputStreams() ;
 892     else
 893       nStreams = fInput ;
 894
 895     Int_t index = 0 ;
 896     for (index = 0 ; index < nStreams ; index++) {
 897       TString tempo(fEventNames[index]) ;
 898       tempo += index ;
 899       printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ;
 900     }
 901
 902  //   printf("\nWith following parameters:\n") ;
 903  //   printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ;
 904  //   printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;
 905  //   printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ;
 906  //   printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;
 907     printf(" ---------------------------------------------------\n") ;
 908   }
 909   else
 910     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
 911
 912 }
 913
 914 //__________________________________________________________________
 915  void AliPHOSDigitizer::PrintDigits(Option_t * option)
 916 {
 917   // Print a table of digits
 918
 919
 920   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 921   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 922   TClonesArray * digits = phosLoader->Digits() ;
 923   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
 924
 925   AliInfo(Form("%d", digits->GetEntriesFast())) ;
 926   printf("\nevent %d", gAlice->GetEvNumber()) ;
 927   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;
 928
 929
 930   if(strstr(option,"all")||strstr(option,"EMC")){
 931     //loop over digits
 932     AliPHOSDigit * digit;
 933     printf("\nEMC digits (with primaries):\n")  ;
 934     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 935     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
 936     Int_t index ;
 937     for (index = 0 ; (index < digits->GetEntriesFast()) &&
 938            (static_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
 939       digit = (AliPHOSDigit * )  digits->At(index) ;
 940       if(digit->GetNprimary() == 0)
 941         continue;
 942 //       printf("%6d  %8d    %6.5e %4d      %2d :",
 943 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
 944       printf("%6d  %.4f    %6.5e %4d      %2d :",
 945               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 946       Int_t iprimary;
 947       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 948         printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 949       }
 950       printf("\n") ;
 951     }
 952   }
 953
 954   if(strstr(option,"all")||strstr(option,"CPV")){
 955
 956     //loop over CPV digits
 957     AliPHOSDigit * digit;
 958     printf("\nCPV digits (with primaries):\n")  ;
 959     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 960     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
 961     Int_t index ;
 962     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
 963       digit = (AliPHOSDigit * )  digits->At(index) ;
 964       if(digit->GetId() > maxEmc){
 965         printf("%6d  %8d    %4d      %2d :",
 966                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
 967         Int_t iprimary;
 968         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 969           printf("%d ",digit->GetPrimary(iprimary+1) ) ;
 970         }
 971         printf("\n") ;
 972       }
 973     }
 974   }
 975
 976 }
 977
 978 //__________________________________________________________________
 979 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
 980 {  // Calculates the time signal generated by noise
 981   //PH  Info("TimeOfNoise", "Change me") ;
 982   return gRandom->Rndm() * 1.28E-5;
 983 }
 984
 985 //__________________________________________________________________
 986 void AliPHOSDigitizer::Unload()
 987 {
 988
 989   Int_t i ;
 990   for(i = 1 ; i < fInput ; i++){
 991     TString tempo(fEventNames[i]) ;
 992     tempo += i ;
 993     AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ;
 994     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 995     phosLoader->UnloadSDigits() ;
 996   }
 997
 998   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
 999   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1000   phosLoader->UnloadDigits() ;
1001 }
1002
1003 //____________________________________________________________________________
1004 void AliPHOSDigitizer::WriteDigits()
1005 {
1006
1007   // Makes TreeD in the output file.
1008   // Check if branch already exists:
1009   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1010   //      already existing branches.
1011   //   else creates branch with Digits, named "PHOS", title "...",
1012   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
1013   //      and names of files, from which digits are made.
1014
1015   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
1016   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1017
1018   const TClonesArray * digits = phosLoader->Digits() ;
1019   TTree * treeD = phosLoader->TreeD();
1020   if(!treeD){
1021     phosLoader->MakeTree("D");
1022     treeD = phosLoader->TreeD();
1023   }
1024
1025   // -- create Digits branch
1026   Int_t bufferSize = 32000 ;
1027   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
1028   digitsBranch->SetTitle(fEventFolderName);
1029   digitsBranch->Fill() ;
1030
1031   phosLoader->WriteDigits("OVERWRITE");
1032
1033   Unload() ;
1034
1035 }