EMCAL/AliEMCALDigitizer.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 //_________________________________________________________________________
  19 //
  20 //////////////////////////////////////////////////////////////////////////////
  21 // Class performs digitization of Summable digits from simulated data
  22 //
  23 // In addition it performs mixing/embedding of summable digits from different events.
  24 //
  25 // For each event 3 branches are created in TreeD:
  26 //   "EMCAL" - list of digits
  27 //   "EMCALTRG" - list of trigger digits
  28 //   "AliEMCALDigitizer" - AliEMCALDigitizer with all parameters used in digitization
  29 //
  30 //
  31 ////////////////////////////////////////////////////////////////////////////////////
  32 //
  33 //*-- Author: Sahal Yacoob (LBL)
  34 // based on : AliEMCALDigitizer
  35 // Modif:
  36 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
  37 //                           of new IO (a la PHOS)
  38 //  November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry
  39 //  July 2011 GCB: Digitizer modified to accomodate embedding.
  40 //                 Time calibration added. Decalibration possibility of energy and time added
  41 //_________________________________________________________________________________
  42
  43 // --- ROOT system ---
  44 #include <TROOT.h>
  45 #include <TTree.h>
  46 #include <TSystem.h>
  47 #include <TBenchmark.h>
  48 #include <TBrowser.h>
  49 #include <TObjectTable.h>
  50 #include <TRandom.h>
  51 #include <TF1.h>
  52 #include <cassert>
  53
  54 // --- AliRoot header files ---
  55 #include "AliLog.h"
  56 #include "AliRun.h"
  57 #include "AliDigitizationInput.h"
  58 #include "AliRunLoader.h"
  59 #include "AliCDBManager.h"
  60 #include "AliCDBEntry.h"
  61 #include "AliEMCALDigit.h"
  62 #include "AliEMCAL.h"
  63 #include "AliEMCALLoader.h"
  64 #include "AliEMCALDigitizer.h"
  65 #include "AliEMCALSDigitizer.h"
  66 #include "AliEMCALGeometry.h"
  67 #include "AliEMCALCalibData.h"
  68 #include "AliEMCALSimParam.h"
  69 #include "AliEMCALRawDigit.h"
  70 #include "AliCaloCalibPedestal.h"
  71
  72   namespace
  73   {
  74     Double_t HeavisideTheta(Double_t x)
  75     {
  76       Double_t signal = 0.;
  77
  78       if (x > 0.) signal = 1.;
  79
  80       return signal;
  81     }
  82
  83     Double_t AnalogFastORFunction(Double_t *x, Double_t *par)
  84     {
  85       Double_t v0 = par[0];
  86       Double_t t0 = par[1];
  87       Double_t tr = par[2];
  88
  89       Double_t R1 = 1000.;
  90       Double_t C1 = 33e-12;
  91       Double_t R2 = 1800;
  92       Double_t C2 = 22e-12;
  93
  94       Double_t t  =   x[0];
  95
  96       return (((0.8*(-((TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-t + t0)/(C1*R1))*
  97                         TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2)) +
  98                      C1*C2*R1*R2*(1 - (C2*TMath::Power(TMath::E(),(-t + t0)/(C2*R2))*R2)/(-(C1*R1) + C2*R2)))*v0*
  99                 HeavisideTheta(t - t0))/tr
 100                - (0.8*(C1*C2*R1*R2 -
 101                        (TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C1*R1))*
 102                         TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2) +
 103                        (C1*TMath::Power(C2,2)*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C2*R2))*
 104                         R1*TMath::Power(R2,2))/(C1*R1 - C2*R2))*v0*
 105                   HeavisideTheta(t - t0 - 1.25*tr))/tr)/(C2*R1));
 106     }
 107   }
 108
 109 ClassImp(AliEMCALDigitizer)
 110
 111
 112 //____________________________________________________________________________
 113 AliEMCALDigitizer::AliEMCALDigitizer()
 114   : AliDigitizer("",""),
 115     fDefaultInit(kTRUE),
 116     fDigitsInRun(0),
 117     fInit(0),
 118     fInput(0),
 119     fInputFileNames(0x0),
 120     fEventNames(0x0),
 121     fDigitThreshold(0),
 122     fMeanPhotonElectron(0),
 123     fGainFluctuations(0),
 124     fPinNoise(0),
 125     fTimeNoise(0),
 126     fTimeDelay(0),
 127     fTimeResolutionPar0(0),
 128     fTimeResolutionPar1(0),
 129     fADCchannelEC(0),
 130     fADCpedestalEC(0),
 131     fADCchannelECDecal(0),
 132     fTimeChannel(0),
 133     fTimeChannelDecal(0),
 134     fNADCEC(0),
 135     fEventFolderName(""),
 136     fFirstEvent(0),
 137     fLastEvent(0),
 138     fCalibData(0x0),
 139     fSDigitizer(0x0)
 140 {
 141   // ctor
 142   InitParameters() ;
 143   fDigInput = 0 ;                     // We work in the standalone mode
 144 }
 145
 146 //____________________________________________________________________________
 147 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
 148   : AliDigitizer("EMCALDigitizer", alirunFileName),
 149     fDefaultInit(kFALSE),
 150     fDigitsInRun(0),
 151     fInit(0),
 152     fInput(0),
 153     fInputFileNames(0),
 154     fEventNames(0),
 155     fDigitThreshold(0),
 156     fMeanPhotonElectron(0),
 157     fGainFluctuations(0),
 158     fPinNoise(0),
 159     fTimeNoise(0),
 160     fTimeDelay(0),
 161     fTimeResolutionPar0(0),
 162     fTimeResolutionPar1(0),
 163     fADCchannelEC(0),
 164     fADCpedestalEC(0),
 165     fADCchannelECDecal(0),
 166     fTimeChannel(0),
 167     fTimeChannelDecal(0),
 168     fNADCEC(0),
 169     fEventFolderName(eventFolderName),
 170     fFirstEvent(0),
 171     fLastEvent(0),
 172     fCalibData(0x0),
 173     fSDigitizer(0x0)
 174 {
 175   // ctor
 176   InitParameters() ;
 177   Init() ;
 178   fDigInput = 0 ;                     // We work in the standalone mode
 179 }
 180
 181 //____________________________________________________________________________
 182 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
 183   : AliDigitizer(d.GetName(),d.GetTitle()),
 184     fDefaultInit(d.fDefaultInit),
 185     fDigitsInRun(d.fDigitsInRun),
 186     fInit(d.fInit),
 187     fInput(d.fInput),
 188     fInputFileNames(d.fInputFileNames),
 189     fEventNames(d.fEventNames),
 190     fDigitThreshold(d.fDigitThreshold),
 191     fMeanPhotonElectron(d.fMeanPhotonElectron),
 192     fGainFluctuations(d.fGainFluctuations),
 193     fPinNoise(d.fPinNoise),
 194     fTimeNoise(d.fTimeNoise),
 195     fTimeDelay(d.fTimeDelay),
 196     fTimeResolutionPar0(d.fTimeResolutionPar0),
 197     fTimeResolutionPar1(d.fTimeResolutionPar1),
 198     fADCchannelEC(d.fADCchannelEC),
 199     fADCpedestalEC(d.fADCpedestalEC),
 200     fADCchannelECDecal(d.fADCchannelECDecal),
 201     fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
 202     fNADCEC(d.fNADCEC),
 203     fEventFolderName(d.fEventFolderName),
 204     fFirstEvent(d.fFirstEvent),
 205     fLastEvent(d.fLastEvent),
 206     fCalibData(d.fCalibData),
 207     fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
 208 {
 209   // copyy ctor
 210 }
 211
 212 //____________________________________________________________________________
 213 AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
 214   : AliDigitizer(rd,"EMCALDigitizer"),
 215     fDefaultInit(kFALSE),
 216     fDigitsInRun(0),
 217     fInit(0),
 218     fInput(0),
 219     fInputFileNames(0),
 220     fEventNames(0),
 221     fDigitThreshold(0),
 222     fMeanPhotonElectron(0),
 223     fGainFluctuations(0),
 224     fPinNoise(0.),
 225     fTimeNoise(0.),
 226     fTimeDelay(0.),
 227     fTimeResolutionPar0(0.),
 228     fTimeResolutionPar1(0.),
 229     fADCchannelEC(0),
 230     fADCpedestalEC(0),
 231     fADCchannelECDecal(0),
 232     fTimeChannel(0),
 233     fTimeChannelDecal(0),
 234     fNADCEC(0),
 235     fEventFolderName(0),
 236     fFirstEvent(0),
 237     fLastEvent(0),
 238     fCalibData(0x0),
 239     fSDigitizer(0x0)
 240 {
 241   // ctor Init() is called by RunDigitizer
 242   fDigInput = rd ;
 243   fEventFolderName = fDigInput->GetInputFolderName(0) ;
 244   SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
 245   InitParameters() ;
 246 }
 247
 248 //____________________________________________________________________________
 249   AliEMCALDigitizer::~AliEMCALDigitizer()
 250 {
 251   //dtor
 252   delete [] fInputFileNames ;
 253   delete [] fEventNames ;
 254   if (fSDigitizer) delete fSDigitizer;
 255 }
 256
 257 //____________________________________________________________________________
 258 void AliEMCALDigitizer::Digitize(Int_t event)
 259 {
 260   // Makes the digitization of the collected summable digits
 261   // for this it first creates the array of all EMCAL modules
 262   // filled with noise and after that adds contributions from
 263   // SDigits. This design helps to avoid scanning over the
 264   // list of digits to add  contribution of any new SDigit.
 265   //
 266   // JLK 26-Jun-2008
 267   // Note that SDigit energy info is stored as an amplitude, so we
 268   // must call the Calibrate() method of the SDigitizer to convert it
 269   // back to an energy in GeV before adding it to the Digit
 270   //
 271
 272   AliRunLoader *rl = AliRunLoader::Instance();
 273   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 274
 275   if(!emcalLoader)
 276   {
 277     AliFatal("EMCALLoader is NULL!") ;
 278     return; // not needed but in case coverity complains ...
 279   }
 280
 281   Int_t readEvent = event ;
 282   if (fDigInput) // fDigInput is data member from AliDigitizer
 283     readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
 284   AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
 285                   readEvent, GetTitle(), fEventFolderName.Data())) ;
 286   rl->GetEvent(readEvent);
 287
 288   TClonesArray * digits = emcalLoader->Digits() ;
 289   digits->Delete() ;  //correct way to clear array when memory is
 290   //allocated by objects stored in array
 291
 292   // Load Geometry
 293   if (!rl->GetAliRun())
 294   {
 295     AliFatal("Could not get AliRun from runLoader");
 296     return; // not needed but in case coverity complains ...
 297   }
 298
 299   AliEMCAL * emcal  = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
 300   AliEMCALGeometry *geom = emcal->GetGeometry();
 301   static int nEMC = geom->GetNCells();//max number of digits possible
 302   AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
 303
 304   digits->Expand(nEMC) ;
 305
 306   // RS create a digitizer on the fly
 307   if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
 308   fSDigitizer->SetEventRange(0, -1) ;
 309
 310   //-------------------------------------------------------
 311   //take all the inputs to add together and load the SDigits
 312   TObjArray * sdigArray = new TObjArray(fInput) ;
 313   sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
 314
 315   Int_t i ;
 316   Int_t embed = kFALSE;
 317   for(i = 1 ; i < fInput ; i++)
 318   {
 319     TString tempo(fEventNames[i]) ;
 320     tempo += i ;
 321
 322     AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ;
 323
 324     if (!rl2)
 325       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
 326
 327     if(!rl2)
 328     {
 329       AliFatal("Run Loader of second event not available!");
 330       return; // not needed but in case coverity complains ...
 331     }
 332
 333     if (fDigInput)
 334       readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
 335
 336     Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
 337
 338     rl2->LoadSDigits();
 339     //     rl2->LoadDigits();
 340     rl2->GetEvent(readEvent);
 341
 342     if(!rl2->GetDetectorLoader("EMCAL"))
 343     {
 344       AliFatal("Loader of second input not found");
 345       return; // not needed but in case coverity complains ...
 346     }
 347
 348     AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
 349
 350     if(!emcalLoader2)
 351     {
 352       AliFatal("EMCAL Loader of second event not available!");
 353       return; // not needed but in case coverity complains ...
 354     }
 355
 356     //Merge 2 simulated sdigits
 357     if(!emcalLoader2->SDigits()) continue;
 358
 359     TClonesArray* sdigits2 = emcalLoader2->SDigits();
 360     sdigArray->AddAt(sdigits2, i) ;
 361
 362     // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
 363     // do not smear energy of embedded, do not add noise to any sdigits
 364     if( sdigits2->GetEntriesFast() <= 0 ) continue;
 365
 366     //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
 367     AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
 368     if( digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded ) embed = kTRUE;
 369
 370   }// input loop
 371
 372   //--------------------------------
 373   //Find the first tower with signal
 374   Int_t nextSig = nEMC + 1 ;
 375   TClonesArray * sdigits ;
 376   for(i = 0 ; i < fInput ; i++)
 377   {
 378     if(i > 0 && embed) continue;
 379
 380     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 381     if(!sdigits)
 382     {
 383       AliDebug(1,"Null sdigit pointer");
 384       continue;
 385     }
 386
 387     if (sdigits->GetEntriesFast() <= 0 )
 388     {
 389       AliDebug(1, "No sdigits entries");
 390       continue;
 391     }
 392
 393     AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
 394     if(!sd)
 395     {
 396       AliDebug(1, "NULL sdigit pointer");
 397       continue;
 398     }
 399
 400     Int_t curNext = sd->GetId() ;
 401     if(curNext < nextSig)
 402     nextSig = curNext ;
 403     AliDebug(1,Form("input %i : #sdigits %i \n",i, sdigits->GetEntriesFast()));
 404
 405   }// input loop
 406
 407   AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
 408
 409   TArrayI index(fInput) ;
 410   index.Reset() ;  //Set all indexes to zero
 411
 412   AliEMCALDigit * digit ;
 413   AliEMCALDigit * curSDigit ;
 414
 415   //---------------------------------------------
 416   //Put Noise contribution, smear time and energy
 417   Float_t timeResolution = 0;
 418   Int_t absID = -1 ;
 419   for(absID = 0; absID < nEMC; absID++)
 420   { // Nov 30, 2006 by PAI; was from 1 to nEMC
 421
 422     if (IsDead(absID)) continue; // Don't instantiate dead digits
 423
 424     Float_t energy = 0 ;
 425
 426     // amplitude set to zero, noise will be added later
 427     Float_t noiseTime = 0.;
 428     if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
 429     new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
 430     //look if we have to add signal?
 431     digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
 432
 433     if (!digit)
 434     {
 435       AliDebug(1,"Digit pointer is null");
 436       continue;
 437     }
 438
 439     if(absID==nextSig)
 440     {
 441       // Calculate time as time of the largest digit
 442       Float_t time = digit->GetTime() ;
 443       Float_t aTime= digit->GetAmplitude() ;
 444
 445       // loop over input
 446       Int_t nInputs = fInput;
 447       if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
 448       for(i = 0; i< nInputs ; i++)
 449       {
 450         //loop over (possible) merge sources
 451         TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
 452         if(sdtclarr)
 453         {
 454           Int_t sDigitEntries = sdtclarr->GetEntriesFast();
 455
 456           if(sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
 457           else                          curSDigit = 0 ;
 458
 459           //May be several digits will contribute from the same input
 460           while(curSDigit && (curSDigit->GetId() == absID))
 461           {
 462             //Shift primary to separate primaries belonging different inputs
 463             Int_t primaryoffset = i ;
 464             if(fDigInput) primaryoffset = fDigInput->GetMask(i) ;
 465             curSDigit->ShiftPrimary(primaryoffset) ;
 466
 467             if(curSDigit->GetAmplitude()>aTime)
 468             {
 469               aTime = curSDigit->GetAmplitude();
 470               time  = curSDigit->GetTime();
 471             }
 472
 473             *digit = *digit + *curSDigit ;  //adds amplitudes of each digit
 474
 475             index[i]++ ;
 476
 477             if( sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
 478             else                           curSDigit = 0 ;
 479           }// while
 480         }// source exists
 481       }// loop over merging sources
 482
 483       //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
 484       energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
 485
 486       // add fluctuations for photo-electron creation
 487       // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
 488       Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
 489       energy       *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
 490
 491       //calculate and set time
 492       digit->SetTime(time) ;
 493
 494       //Find next signal module
 495       nextSig = nEMC + 1 ;
 496       for(i = 0 ; i < nInputs ; i++)
 497       {
 498         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 499
 500         if(sdigits)
 501         {
 502           Int_t curNext = nextSig ;
 503           if(sdigits->GetEntriesFast() > index[i])
 504           {
 505             AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
 506             if(tmpdigit)
 507             {
 508               curNext = tmpdigit->GetId() ;
 509             }
 510           }
 511           if(curNext < nextSig) nextSig = curNext ;
 512         }// sdigits exist
 513       } // input loop
 514
 515     }//absID==nextSig
 516
 517     // add the noise now, no need for embedded events with real data
 518     if(!embed)
 519       energy += gRandom->Gaus(0., fPinNoise) ;
 520
 521     // JLK 26-June-2008
 522     //Now digitize the energy according to the fSDigitizer method,
 523     //which merely converts the energy to an integer.  Later we will
 524     //check that the stored value matches our allowed dynamic ranges
 525     digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
 526
 527     //Set time resolution with final energy
 528     timeResolution = GetTimeResolution(energy);
 529     digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
 530     AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
 531                      absID, energy, nextSig));
 532     //Add delay to time
 533     digit->SetTime(digit->GetTime()+fTimeDelay) ;
 534
 535   } // for(absID = 0; absID < nEMC; absID++)
 536
 537   //---------------------------------------------------------
 538   //Embed simulated amplitude (and time?) to real data digits
 539   if(embed)
 540   {
 541     for(Int_t input = 1; input<fInput; input++)
 542     {
 543       TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
 544       if(!realDigits)
 545       {
 546         AliDebug(1,"No sdigits to merge\n");
 547         continue;
 548       }
 549
 550       for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
 551       {
 552         AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
 553         if(digit2)
 554         {
 555           digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
 556           if(digit)
 557           {
 558             // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
 559             // and multiply to get a big int.
 560             Float_t time2 = digit2->GetTime();
 561             Float_t e2    = digit2->GetAmplitude();
 562             CalibrateADCTime(e2,time2,digit2->GetId());
 563             Float_t amp2 = fSDigitizer->Digitize(e2);
 564
 565             Float_t e0    = digit ->GetAmplitude();
 566             if(e0>0.01)
 567             AliDebug(1,Form("digit 1: Abs ID %d, amp %f, type %d, time %e; digit2: Abs ID %d, amp %f, type %d, time %e\n",
 568                             digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
 569                             digit2->GetId(),amp2,                   digit2->GetType(), time2           ));
 570
 571             // Sum amplitudes, change digit amplitude
 572             digit->SetAmplitude( digit->GetAmplitude() + amp2);
 573             // Rough assumption, assign time of the largest of the 2 digits,
 574             // data or signal to the final digit.
 575             if(amp2 > digit->GetAmplitude())  digit->SetTime(time2);
 576             //Mark the new digit as embedded
 577             digit->SetType(AliEMCALDigit::kEmbedded);
 578
 579             if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
 580             AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
 581                             digit->GetId(), digit->GetAmplitude(), digit->GetType()));
 582           }//digit2
 583         }//digit2
 584       }//loop digit 2
 585     }//input loop
 586   }//embed
 587
 588   //JLK is it better to call Clear() here?
 589   delete sdigArray ; //We should not delete its contents
 590
 591   //------------------------------
 592   //remove digits below thresholds
 593   // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
 594   // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
 595   // before merge in the same loop real data digits if available
 596   Float_t energy = 0;
 597   Float_t time   = 0;
 598   for(i = 0 ; i < nEMC ; i++)
 599   {
 600     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
 601     if(digit)
 602     {
 603       //Then get the energy in GeV units.
 604       energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
 605       //Then digitize using the calibration constants of the ocdb
 606       Float_t ampADC = energy;
 607       DigitizeEnergyTime(ampADC, time, digit->GetId())  ;
 608       if(ampADC < fDigitThreshold)
 609       digits->RemoveAt(i) ;
 610
 611     }// digit exists
 612   } // digit loop
 613
 614   digits->Compress() ;
 615
 616   Int_t ndigits = digits->GetEntriesFast() ;
 617
 618   //---------------------------------------------------------------
 619   //JLK 26-June-2008
 620   //After we have done the summing and digitizing to create the
 621   //digits, now we want to calibrate the resulting amplitude to match
 622   //the dynamic range of our real data.
 623   for (i = 0 ; i < ndigits ; i++)
 624   {
 625     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
 626     if(digit)
 627     {
 628       digit->SetIndexInList(i) ;
 629       energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
 630       time   = digit->GetTime();
 631       Float_t ampADC = energy;
 632       DigitizeEnergyTime(ampADC, time, digit->GetId());
 633       digit->SetAmplitude(ampADC) ;
 634       digit->SetTime(time);
 635       // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
 636     }// digit exists
 637   }//Digit loop
 638
 639 }
 640
 641 //_____________________________________________________________________
 642 void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
 643 {
 644   // JLK 26-June-2008
 645   // Returns digitized value of the energy in a cell absId
 646   // using the calibration constants stored in the OCDB
 647   // or default values if no CalibData object is found.
 648   // This effectively converts everything to match the dynamic range
 649   // of the real data we will collect
 650   //
 651   // Load Geometry
 652   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 653
 654   if (geom==0)
 655   {
 656     AliFatal("Did not get geometry from EMCALLoader");
 657     return;
 658   }
 659
 660   Int_t iSupMod = -1;
 661   Int_t nModule = -1;
 662   Int_t nIphi   = -1;
 663   Int_t nIeta   = -1;
 664   Int_t iphi    = -1;
 665   Int_t ieta    = -1;
 666   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 667   if(!bCell)
 668   Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
 669   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 670
 671   if(fCalibData)
 672   {
 673     fADCpedestalEC     = fCalibData->GetADCpedestal     (iSupMod,ieta,iphi);
 674     fADCchannelEC      = fCalibData->GetADCchannel      (iSupMod,ieta,iphi);
 675     fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
 676     fTimeChannel       = fCalibData->GetTimeChannel     (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
 677     fTimeChannelDecal  = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
 678   }
 679
 680   //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
 681   energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal   ;
 682   time  += fTimeChannel-fTimeChannelDecal;
 683
 684   if(energy > fNADCEC )
 685   energy =  fNADCEC ;
 686 }
 687
 688 //_____________________________________________________________________
 689 void AliEMCALDigitizer::Decalibrate(AliEMCALDigit *digit)
 690 {
 691   // Decalibrate, used in Trigger digits
 692   // Load Geometry
 693   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 694
 695   if (!geom)
 696   {
 697     AliFatal("Did not get geometry from EMCALLoader");
 698     return;
 699   }
 700
 701   Int_t absId   = digit->GetId();
 702   Int_t iSupMod = -1;
 703   Int_t nModule = -1;
 704   Int_t nIphi   = -1;
 705   Int_t nIeta   = -1;
 706   Int_t iphi    = -1;
 707   Int_t ieta    = -1;
 708
 709   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 710
 711   if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
 712   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 713
 714   if (fCalibData)
 715   {
 716     fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
 717     Float_t factor = 1./(fADCchannelEC/0.0162);
 718
 719     *digit = *digit * factor;
 720   }
 721 }
 722
 723 //_____________________________________________________________________
 724 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
 725 {
 726   // Returns the energy in a cell absId with a given adc value
 727   // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
 728   // Used in case of embedding, transform ADC counts from real event into energy
 729   // so that we can add the energy of the simulated sdigits which are in energy
 730   // units.
 731   // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
 732   // or with time out of window
 733
 734   // Load Geometry
 735   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 736
 737   if (!geom)
 738   {
 739     AliFatal("Did not get geometry from EMCALLoader");
 740     return;
 741   }
 742
 743   Int_t iSupMod = -1;
 744   Int_t nModule = -1;
 745   Int_t nIphi   = -1;
 746   Int_t nIeta   = -1;
 747   Int_t iphi    = -1;
 748   Int_t ieta    = -1;
 749   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 750   if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
 751   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 752
 753   if(fCalibData)
 754   {
 755     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 756     fADCchannelEC  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 757     fTimeChannel   = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
 758   }
 759
 760   adc   = adc * fADCchannelEC - fADCpedestalEC;
 761   time -= fTimeChannel;
 762
 763 }
 764
 765
 766 //____________________________________________________________________________
 767 void AliEMCALDigitizer::Digitize(Option_t *option)
 768 {
 769   // Steering method to process digitization for events
 770   // in the range from fFirstEvent to fLastEvent.
 771   // This range is optionally set by SetEventRange().
 772   // if fLastEvent=-1, then process events until the end.
 773   // by default fLastEvent = fFirstEvent (process only one event)
 774
 775   if (!fInit) { // to prevent overwrite existing file
 776     Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
 777     return ;
 778   }
 779
 780   if (strstr(option,"print")) {
 781
 782     Print();
 783     return ;
 784   }
 785
 786   if(strstr(option,"tim"))
 787     gBenchmark->Start("EMCALDigitizer");
 788
 789   AliRunLoader *rl = AliRunLoader::Instance();
 790   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 791   Int_t nEvents = 0;
 792   if(!emcalLoader){
 793     AliFatal("Did not get the  Loader");
 794   }
 795   else{
 796
 797     if (fLastEvent == -1)  {
 798       fLastEvent = rl->GetNumberOfEvents() - 1 ;
 799     }
 800     else if (fDigInput)
 801       fLastEvent = fFirstEvent ; // what is this ??
 802
 803     nEvents = fLastEvent - fFirstEvent + 1;
 804     Int_t ievent  = -1;
 805
 806     AliEMCALGeometry *geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
 807     const Int_t nTRU = geom->GetNTotalTRU();
 808     TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit",    nTRU*96);
 809     TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
 810
 811     rl->LoadSDigits("EMCAL");
 812     for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 813
 814       rl->GetEvent(ievent);
 815
 816       Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 817
 818       WriteDigits() ;
 819
 820       //Trigger Digits
 821       //-------------------------------------
 822
 823
 824       Digits2FastOR(digitsTMP, digitsTRG);
 825
 826       WriteDigits(digitsTRG);
 827
 828       (emcalLoader->TreeD())->Fill();
 829
 830       emcalLoader->WriteDigits(   "OVERWRITE");
 831
 832       Unload();
 833
 834       digitsTRG  ->Delete();
 835       digitsTMP  ->Delete();
 836
 837       //-------------------------------------
 838
 839       if(strstr(option,"deb"))
 840         PrintDigits(option);
 841       if(strstr(option,"table")) gObjectTable->Print();
 842
 843       //increment the total number of Digits per run
 844       fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
 845     }//loop
 846
 847   }//loader exists
 848
 849   if(strstr(option,"tim")){
 850     gBenchmark->Stop("EMCALDigitizer");
 851     Float_t cputime   = gBenchmark->GetCpuTime("EMCALDigitizer");
 852     Float_t avcputime = cputime;
 853     if(nEvents==0) avcputime = 0 ;
 854     AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
 855   }
 856 }
 857
 858 //__________________________________________________________________
 859 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
 860 {
 861   // Assign a smeared time to the digit, from observed distribution in LED system (?).
 862   // From F. Blanco
 863   Float_t res = -1;
 864   if (energy > 0)
 865   {
 866     res = TMath::Sqrt(fTimeResolutionPar0 +
 867                       fTimeResolutionPar1/(energy*energy) );
 868   }
 869   // parametrization above is for ns. Convert to seconds:
 870   return res*1e-9;
 871 }
 872
 873 //____________________________________________________________________________
 874 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
 875 {
 876   // FEE digits afterburner to produce TRG digits
 877   // we are only interested in the FEE digit deposited energy
 878   // to be converted later into a voltage value
 879
 880   // push the FEE digit to its associated FastOR (numbered from 0:95)
 881   // TRU is in charge of summing module digits
 882
 883   AliRunLoader *runLoader = AliRunLoader::Instance();
 884
 885   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
 886   if (!emcalLoader) AliFatal("Did not get the  Loader");
 887
 888   const  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 889
 890   // build FOR from simulated digits
 891   // and xfer to the corresponding TRU input (mapping)
 892
 893   TClonesArray* sdigits = emcalLoader->SDigits();
 894
 895   TClonesArray *digits = (TClonesArray*)sdigits->Clone();
 896
 897   AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
 898
 899   TIter NextDigit(digits);
 900   while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
 901   {
 902     if (IsDead(digit)) continue;
 903
 904     Decalibrate(digit);
 905
 906     Int_t id = digit->GetId();
 907
 908     Int_t trgid;
 909     if (geom && geom->GetFastORIndexFromCellIndex(id , trgid))
 910     {
 911       AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
 912
 913       AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
 914
 915       if (!d)
 916       {
 917         new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
 918         d = (AliEMCALDigit*)digitsTMP->At(trgid);
 919         d->SetId(trgid);
 920       }
 921       else
 922       {
 923         *d = *d + *digit;
 924       }
 925     }
 926   }
 927
 928   if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
 929
 930   Int_t     nSamples = geom->GetNTotalTRU();
 931   Int_t *timeSamples = new Int_t[nSamples];
 932
 933   NextDigit = TIter(digitsTMP);
 934   while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
 935   {
 936     if (digit)
 937     {
 938       Int_t     id = digit->GetId();
 939       Float_t time = 50.e-9;
 940
 941       Double_t depositedEnergy = 0.;
 942       for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
 943
 944       if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
 945
 946       // FIXME: Check digit time!
 947       if (depositedEnergy) {
 948         depositedEnergy += gRandom->Gaus(0., .08);
 949         DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
 950
 951         for (Int_t j=0;j<nSamples;j++) {
 952           if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
 953           timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
 954         }
 955
 956         new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
 957
 958         if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
 959       }
 960     }
 961   }
 962
 963   delete [] timeSamples;
 964
 965   if (digits && digits->GetEntriesFast()) digits->Delete();
 966 }
 967
 968 //____________________________________________________________________________
 969 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
 970 {
 971   // parameters:
 972   // id: 0..95
 973   const Int_t    reso = 12;      // 11-bit resolution ADC
 974   const Double_t vFSR = 2.;      // Full scale input voltage range 2V (p-p)
 975 //const Double_t dNe  = 125;     // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
 976   const Double_t dNe  = 125/1.3; // F-ALTRO max V. FEE: factor 4
 977   const Double_t vA   = .136e-6; // CSP output range: 0.136uV/e-
 978   const Double_t rise = 50e-9;   // rise time (10-90%) of the FastOR signal before shaping
 979
 980   const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
 981
 982   Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
 983
 984   TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
 985   signalF.SetParameter( 0,   vV );
 986   signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
 987   signalF.SetParameter( 2, rise );
 988
 989   for (Int_t iTime=0; iTime<nSamples; iTime++)
 990   {
 991     // FIXME: add noise (probably not simply Gaussian) according to DA measurements
 992     // probably plan an access to OCDB
 993     Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
 994     if (TMath::Abs(sig) > vFSR/2.) {
 995       AliError("Signal overflow!");
 996       timeSamples[iTime] = (1 << reso) - 1;
 997     } else {
 998       AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
 999       timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
1000     }
1001   }
1002 }
1003
1004 //____________________________________________________________________________
1005 Bool_t AliEMCALDigitizer::Init()
1006 {
1007   // Makes all memory allocations
1008   fInit = kTRUE ;
1009   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1010
1011   if ( emcalLoader == 0 ) {
1012     Fatal("Init", "Could not obtain the AliEMCALLoader");
1013     return kFALSE;
1014   }
1015
1016   fFirstEvent = 0 ;
1017   fLastEvent = fFirstEvent ;
1018
1019   if(fDigInput)
1020     fInput = fDigInput->GetNinputs() ;
1021   else
1022     fInput           = 1 ;
1023
1024   fInputFileNames  = new TString[fInput] ;
1025   fEventNames      = new TString[fInput] ;
1026   fInputFileNames[0] = GetTitle() ;
1027   fEventNames[0]     = fEventFolderName.Data() ;
1028   Int_t index ;
1029   for (index = 1 ; index < fInput ; index++) {
1030     fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1031     TString tempo = fDigInput->GetInputFolderName(index) ;
1032     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
1033   }
1034
1035   //Calibration instance
1036   fCalibData = emcalLoader->CalibData();
1037   return fInit ;
1038 }
1039
1040 //____________________________________________________________________________
1041 void AliEMCALDigitizer::InitParameters()
1042 {
1043   // Parameter initialization for digitizer
1044
1045   // Get the parameters from the OCDB via the loader
1046   AliRunLoader *rl = AliRunLoader::Instance();
1047   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1048   AliEMCALSimParam * simParam = 0x0;
1049   if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1050
1051   if(!simParam){
1052     simParam = AliEMCALSimParam::GetInstance();
1053     AliWarning("Simulation Parameters not available in OCDB?");
1054   }
1055
1056   fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400;  // electrons per GeV
1057   fGainFluctuations   = simParam->GetGainFluctuations()   ; // 15.0;
1058
1059   fPinNoise           = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1060   if (fPinNoise < 0.0001 )
1061     Warning("InitParameters", "No noise added\n") ;
1062   fTimeNoise          = simParam->GetTimeNoise(); // 1.28E-5 s
1063   fDigitThreshold     = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1064   fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1065   fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1066   fTimeDelay          = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1067
1068   // These defaults are normally not used.
1069   // Values are read from calibration database instead
1070   fADCchannelEC       = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
1071   fADCpedestalEC      = 0.0 ;   // GeV
1072   fADCchannelECDecal  = 1.0;    // No decalibration by default
1073   fTimeChannel        = 0.0;    // No time calibration by default
1074   fTimeChannelDecal   = 0.0;    // No time decalibration by default
1075
1076   fNADCEC             = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
1077
1078   AliDebug(2,Form("Mean Photon Electron %d, Gain Fluct. %2.1f; Noise: APD %f, Time %f; Digit Threshold %d,Time Resolution Par0 %g Par1 %g,NADCEC %d",
1079                   fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1080
1081 }
1082
1083 //__________________________________________________________________
1084 void AliEMCALDigitizer::Print1(Option_t * option)
1085 { // 19-nov-04 - just for convenience
1086   Print();
1087   PrintDigits(option);
1088 }
1089
1090 //__________________________________________________________________
1091 void AliEMCALDigitizer::Print (Option_t * ) const
1092 {
1093   // Print Digitizer's parameters
1094   printf("Print: \n------------------- %s -------------", GetName() ) ;
1095   if( strcmp(fEventFolderName.Data(), "") != 0 ){
1096     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
1097
1098     Int_t nStreams ;
1099     if (fDigInput)
1100       nStreams =  GetNInputStreams() ;
1101     else
1102       nStreams = fInput ;
1103
1104     AliRunLoader *rl=0;
1105
1106     Int_t index = 0 ;
1107     for (index = 0 ; index < nStreams ; index++) {
1108       TString tempo(fEventNames[index]) ;
1109       tempo += index ;
1110       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1111         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1112       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1113       if(emcalLoader){
1114         TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1115         if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1116           fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
1117         printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1118       }//loader
1119     }
1120
1121     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1122
1123     if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1124     else printf("\nNULL LOADER");
1125
1126     printf("\nWith following parameters:\n") ;
1127     printf("    Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1128     printf("    Threshold  in Tower  (fDigitThreshold) = %d\n", fDigitThreshold)  ;
1129     printf("---------------------------------------------------\n")  ;
1130   }
1131   else
1132     printf("Print: AliEMCALDigitizer not initialized") ;
1133 }
1134
1135 //__________________________________________________________________
1136 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1137 {
1138   //utility method for printing digit information
1139
1140   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1141   if(emcalLoader){
1142     TClonesArray * digits  = emcalLoader->Digits() ;
1143     TClonesArray * sdigits = emcalLoader->SDigits() ;
1144
1145     printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1146     printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1147
1148     if(strstr(option,"all")){
1149       //loop over digits
1150       AliEMCALDigit * digit;
1151       printf("\nEMC digits (with primaries):\n")  ;
1152       printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
1153       Int_t index ;
1154       for (index = 0 ; index < digits->GetEntries() ; index++) {
1155         digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1156         if(digit){
1157           printf("\n%6d  %8f    %6.5e %4d      %2d : ",
1158                  digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1159           Int_t iprimary;
1160           for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1161             printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1162           }
1163         }// digit exists
1164       }// loop
1165     }
1166     printf("\n");
1167   }// loader exists
1168   else printf("NULL LOADER, cannot print\n");
1169 }
1170
1171 //__________________________________________________________________
1172 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1173 {
1174   // Calculates the time signal generated by noise
1175   //printf("Time noise %e\n",fTimeNoise);
1176   return gRandom->Rndm() * fTimeNoise;
1177 }
1178
1179 //__________________________________________________________________
1180 void AliEMCALDigitizer::Unload()
1181 {
1182   // Unloads the SDigits and Digits
1183   AliRunLoader *rl=0;
1184
1185   Int_t i ;
1186   for(i = 1 ; i < fInput ; i++){
1187     TString tempo(fEventNames[i]) ;
1188     tempo += i;
1189     if ((rl = AliRunLoader::GetRunLoader(tempo)))
1190       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1191   }
1192   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1193   if(emcalLoader)emcalLoader->UnloadDigits() ;
1194   else AliFatal("Did not get the loader");
1195 }
1196
1197 //_________________________________________________________________________________________
1198 void AliEMCALDigitizer::WriteDigits()
1199 { // Makes TreeD in the output file.
1200   // Check if branch already exists:
1201   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1202   //      already existing branches.
1203   //   else creates branch with Digits, named "EMCAL", title "...",
1204   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1205   //      and names of files, from which digits are made.
1206
1207   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1208
1209   if(emcalLoader){
1210     const TClonesArray * digits = emcalLoader->Digits() ;
1211     TTree * treeD = emcalLoader->TreeD();
1212     if ( !treeD ) {
1213       emcalLoader->MakeDigitsContainer();
1214       treeD = emcalLoader->TreeD();
1215     }
1216
1217     // -- create Digits branch
1218     Int_t bufferSize = 32000 ;
1219     TBranch * digitsBranch = 0;
1220     if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1221       digitsBranch->SetAddress(&digits);
1222       AliWarning("Digits Branch already exists. Not all digits will be visible");
1223     }
1224     else
1225       treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1226     //digitsBranch->SetTitle(fEventFolderName);
1227
1228     //  treeD->Fill() ;
1229     /*
1230      emcalLoader->WriteDigits("OVERWRITE");
1231      emcalLoader->WriteDigitizer("OVERWRITE");
1232
1233      Unload() ;
1234      */
1235
1236   }// loader exists
1237   else AliFatal("Loader not available");
1238 }
1239
1240 //__________________________________________________________________
1241 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1242 { // overloaded method
1243   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1244   if(emcalLoader){
1245
1246     TTree* treeD = emcalLoader->TreeD();
1247     if (!treeD)
1248       {
1249         emcalLoader->MakeDigitsContainer();
1250         treeD = emcalLoader->TreeD();
1251       }
1252
1253     // -- create Digits branch
1254     Int_t bufferSize = 32000;
1255
1256     if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1257       {
1258         triggerBranch->SetAddress(&digits);
1259       }
1260     else
1261       {
1262         treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1263       }
1264
1265     //  treeD->Fill();
1266   }// loader exists
1267   else AliFatal("Loader not available");
1268 }
1269
1270 //__________________________________________________________________
1271 Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1272 {
1273   AliRunLoader *runLoader = AliRunLoader::Instance();
1274   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1275   if (!emcalLoader) AliFatal("Did not get the  Loader");
1276
1277   AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1278   if (!caloPed) {
1279     AliWarning("Could not access pedestal data! No dead channel removal applied");
1280     return kFALSE;
1281   }
1282
1283         // Load Geometry
1284   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1285   if (!geom) AliFatal("Did not get geometry from EMCALLoader");
1286
1287   Int_t absId   = digit->GetId();
1288   Int_t iSupMod = -1;
1289   Int_t nModule = -1;
1290   Int_t nIphi   = -1;
1291   Int_t nIeta   = -1;
1292   Int_t iphi    = -1;
1293   Int_t ieta    = -1;
1294
1295   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1296
1297   if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1298   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1299
1300   Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1301
1302   if (channelStatus == AliCaloCalibPedestal::kDead)
1303     return kTRUE;
1304   else
1305     return kFALSE;
1306 }
1307
1308
1309 //__________________________________________________________________
1310 Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1311 {
1312     AliRunLoader *runLoader = AliRunLoader::Instance();
1313     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1314     if (!emcalLoader) AliFatal("Did not get the  Loader");
1315
1316     AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1317     if (!caloPed) {
1318         AliWarning("Could not access pedestal data! No dead channel removal applied");
1319         return kFALSE;
1320     }
1321
1322         // Load Geometry
1323     const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1324     if (!geom) AliFatal("Did not get geometry from EMCALLoader");
1325
1326     Int_t iSupMod = -1;
1327     Int_t nModule = -1;
1328     Int_t nIphi   = -1;
1329     Int_t nIeta   = -1;
1330     Int_t iphi    = -1;
1331     Int_t ieta    = -1;
1332
1333     Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1334
1335     if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1336     geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1337
1338     Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1339
1340     if (channelStatus == AliCaloCalibPedestal::kDead)
1341         return kTRUE;
1342     else
1343         return kFALSE;
1344 }
1345