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   {
 421     Float_t energy = 0 ;
 422
 423     // amplitude set to zero, noise will be added later
 424     Float_t noiseTime = 0.;
 425     if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
 426     new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
 427     //look if we have to add signal?
 428     digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
 429
 430     if (!digit)
 431     {
 432       AliDebug(1,"Digit pointer is null");
 433       continue;
 434     }
 435
 436     if(absID==nextSig)
 437     {
 438       // Calculate time as time of the largest digit
 439       Float_t time = digit->GetTime() ;
 440       Float_t aTime= digit->GetAmplitude() ;
 441
 442       // loop over input
 443       Int_t nInputs = fInput;
 444       if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
 445       for(i = 0; i< nInputs ; i++)
 446       {
 447         //loop over (possible) merge sources
 448         TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
 449         if(sdtclarr)
 450         {
 451           Int_t sDigitEntries = sdtclarr->GetEntriesFast();
 452
 453           if(sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
 454           else                          curSDigit = 0 ;
 455
 456           //May be several digits will contribute from the same input
 457           while(curSDigit && (curSDigit->GetId() == absID))
 458           {
 459             //Shift primary to separate primaries belonging different inputs
 460             Int_t primaryoffset = i ;
 461             if(fDigInput) primaryoffset = fDigInput->GetMask(i) ;
 462             curSDigit->ShiftPrimary(primaryoffset) ;
 463
 464             if(curSDigit->GetAmplitude()>aTime)
 465             {
 466               aTime = curSDigit->GetAmplitude();
 467               time  = curSDigit->GetTime();
 468             }
 469
 470             *digit = *digit + *curSDigit ;  //adds amplitudes of each digit
 471
 472             index[i]++ ;
 473
 474             if( sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
 475             else                           curSDigit = 0 ;
 476           }// while
 477         }// source exists
 478       }// loop over merging sources
 479
 480       //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
 481       energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
 482
 483       // add fluctuations for photo-electron creation
 484       // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
 485       Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
 486       energy       *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
 487
 488       //calculate and set time
 489       digit->SetTime(time) ;
 490
 491       //Find next signal module
 492       nextSig = nEMC + 1 ;
 493       for(i = 0 ; i < nInputs ; i++)
 494       {
 495         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
 496
 497         if(sdigits)
 498         {
 499           Int_t curNext = nextSig ;
 500           if(sdigits->GetEntriesFast() > index[i])
 501           {
 502             AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
 503             if ( tmpdigit ) curNext = tmpdigit->GetId() ;
 504           }
 505
 506           if(curNext < nextSig) nextSig = curNext ;
 507         }// sdigits exist
 508       } // input loop
 509
 510     }//absID==nextSig
 511
 512     // add the noise now, no need for embedded events with real data
 513     if(!embed)
 514       energy += gRandom->Gaus(0., fPinNoise) ;
 515
 516     // JLK 26-June-2008
 517     //Now digitize the energy according to the fSDigitizer method,
 518     //which merely converts the energy to an integer.  Later we will
 519     //check that the stored value matches our allowed dynamic ranges
 520     digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
 521
 522     //Set time resolution with final energy
 523     timeResolution = GetTimeResolution(energy);
 524     digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
 525     AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
 526                      absID, energy, nextSig));
 527     //Add delay to time
 528     digit->SetTime(digit->GetTime()+fTimeDelay) ;
 529
 530   } // for(absID = 0; absID < nEMC; absID++)
 531
 532   //---------------------------------------------------------
 533   //Embed simulated amplitude (and time?) to real data digits
 534   if ( embed )
 535   {
 536     for(Int_t input = 1; input<fInput; input++)
 537     {
 538       TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
 539       if(!realDigits)
 540       {
 541         AliDebug(1,"No sdigits to merge\n");
 542         continue;
 543       }
 544
 545       for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
 546       {
 547         AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
 548         if ( !digit2 ) continue;
 549
 550         digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
 551         if ( !digit ) continue;
 552
 553         // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
 554         // and multiply to get a big int.
 555         Float_t time2 = digit2->GetTime();
 556         Float_t e2    = digit2->GetAmplitude();
 557         CalibrateADCTime(e2,time2,digit2->GetId());
 558         Float_t amp2  = fSDigitizer->Digitize(e2);
 559
 560         Float_t e0    = digit ->GetAmplitude();
 561         if(e0>0.01)
 562           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",
 563                           digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
 564                           digit2->GetId(),amp2,                   digit2->GetType(), time2           ));
 565
 566         // Sum amplitudes, change digit amplitude
 567         digit->SetAmplitude( digit->GetAmplitude() + amp2);
 568
 569         // Rough assumption, assign time of the largest of the 2 digits,
 570         // data or signal to the final digit.
 571         if(amp2 > digit->GetAmplitude())  digit->SetTime(time2);
 572
 573         //Mark the new digit as embedded
 574         digit->SetType(AliEMCALDigit::kEmbedded);
 575
 576         if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
 577           AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
 578                           digit->GetId(), digit->GetAmplitude(), digit->GetType()));
 579       }//loop digit 2
 580     }//input loop
 581   }//embed
 582
 583   //JLK is it better to call Clear() here?
 584   delete sdigArray ; //We should not delete its contents
 585
 586   //------------------------------
 587   //remove digits below thresholds
 588   // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
 589   // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
 590   // before merge in the same loop real data digits if available
 591   Float_t energy = 0;
 592   Float_t time   = 0;
 593   for(i = 0 ; i < nEMC ; i++)
 594   {
 595     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
 596     if ( !digit ) continue;
 597
 598     //Then get the energy in GeV units.
 599     energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
 600
 601     //Then digitize using the calibration constants of the ocdb
 602     Float_t ampADC = energy;
 603     DigitizeEnergyTime(ampADC, time, digit->GetId())  ;
 604
 605     if(ampADC < fDigitThreshold || IsDead(digit->GetId()))
 606       digits->RemoveAt(i) ;
 607
 608   } // digit loop
 609
 610   digits->Compress() ;
 611
 612   Int_t ndigits = digits->GetEntriesFast() ;
 613
 614   //---------------------------------------------------------------
 615   //JLK 26-June-2008
 616   //After we have done the summing and digitizing to create the
 617   //digits, now we want to calibrate the resulting amplitude to match
 618   //the dynamic range of our real data.
 619   for (i = 0 ; i < ndigits ; i++)
 620   {
 621     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
 622     if( !digit ) continue ;
 623
 624     digit->SetIndexInList(i) ;
 625
 626     time   = digit->GetTime();
 627     digit->SetTime(time);
 628
 629     energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
 630
 631     Float_t ampADC = energy;
 632     DigitizeEnergyTime(ampADC, time, digit->GetId());
 633
 634     digit->SetAmplitude(ampADC) ;
 635     // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
 636   }//Digit loop
 637
 638 }
 639
 640 //_____________________________________________________________________
 641 void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
 642 {
 643   // JLK 26-June-2008
 644   // Returns digitized value of the energy in a cell absId
 645   // using the calibration constants stored in the OCDB
 646   // or default values if no CalibData object is found.
 647   // This effectively converts everything to match the dynamic range
 648   // of the real data we will collect
 649   //
 650   // Load Geometry
 651   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 652
 653   if (geom==0)
 654   {
 655     AliFatal("Did not get geometry from EMCALLoader");
 656     return;
 657   }
 658
 659   Int_t iSupMod = -1;
 660   Int_t nModule = -1;
 661   Int_t nIphi   = -1;
 662   Int_t nIeta   = -1;
 663   Int_t iphi    = -1;
 664   Int_t ieta    = -1;
 665   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 666   if(!bCell)
 667   Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
 668   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 669
 670   if(fCalibData)
 671   {
 672     fADCpedestalEC     = fCalibData->GetADCpedestal     (iSupMod,ieta,iphi);
 673     fADCchannelEC      = fCalibData->GetADCchannel      (iSupMod,ieta,iphi);
 674     fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
 675     fTimeChannel       = fCalibData->GetTimeChannel     (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
 676     fTimeChannelDecal  = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
 677   }
 678
 679   //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
 680   energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal   ;
 681   time  += fTimeChannel-fTimeChannelDecal;
 682
 683   if ( energy > fNADCEC ) energy =  fNADCEC ;
 684 }
 685
 686 //_____________________________________________________________________
 687 void AliEMCALDigitizer::DecalibrateTrigger(AliEMCALDigit *digit)
 688 {
 689   // Decalibrate, used in Trigger digits
 690
 691   if ( !fCalibData ) return ;
 692
 693   // Load Geometry
 694   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 695
 696   if (!geom)
 697   {
 698     AliFatal("Did not get geometry from EMCALLoader");
 699     return;
 700   }
 701
 702   // Recover the digit location in row-column-SM
 703   Int_t absId   = digit->GetId();
 704   Int_t iSupMod = -1;
 705   Int_t nModule = -1;
 706   Int_t nIphi   = -1;
 707   Int_t nIeta   = -1;
 708   Int_t iphi    = -1;
 709   Int_t ieta    = -1;
 710
 711   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 712   if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
 713
 714   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 715
 716   // Now decalibrate
 717   Float_t adcChannel       = fCalibData->GetADCchannel      (iSupMod,ieta,iphi);
 718   Float_t adcChannelOnline = fCalibData->GetADCchannelOnline(iSupMod,ieta,iphi);
 719   //printf("calib param for (SM%d,col%d,row%d): %1.4f - online param: %1.4f \n",iSupMod,ieta,iphi,adcChannel,adcChannelOnline);
 720   Float_t factor = 1./(adcChannel/adcChannelOnline);
 721
 722   *digit = *digit * factor;
 723
 724 }
 725
 726 //_____________________________________________________________________
 727 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
 728 {
 729   // Returns the energy in a cell absId with a given adc value
 730   // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
 731   // Used in case of embedding, transform ADC counts from real event into energy
 732   // so that we can add the energy of the simulated sdigits which are in energy
 733   // units.
 734   // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
 735   // or with time out of window
 736
 737   // Load Geometry
 738   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 739
 740   if (!geom)
 741   {
 742     AliFatal("Did not get geometry from EMCALLoader");
 743     return;
 744   }
 745
 746   Int_t iSupMod = -1;
 747   Int_t nModule = -1;
 748   Int_t nIphi   = -1;
 749   Int_t nIeta   = -1;
 750   Int_t iphi    = -1;
 751   Int_t ieta    = -1;
 752   Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 753   if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
 754   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 755
 756   if(fCalibData)
 757   {
 758     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 759     fADCchannelEC  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 760     fTimeChannel   = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
 761   }
 762
 763   adc   = adc * fADCchannelEC - fADCpedestalEC;
 764   time -= fTimeChannel;
 765
 766 }
 767
 768
 769 //____________________________________________________________________________
 770 void AliEMCALDigitizer::Digitize(Option_t *option)
 771 {
 772   // Steering method to process digitization for events
 773   // in the range from fFirstEvent to fLastEvent.
 774   // This range is optionally set by SetEventRange().
 775   // if fLastEvent=-1, then process events until the end.
 776   // by default fLastEvent = fFirstEvent (process only one event)
 777
 778   if (!fInit)
 779   { // to prevent overwrite existing file
 780     Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
 781     return ;
 782   }
 783
 784   if (strstr(option,"print"))
 785   {
 786     Print();
 787     return ;
 788   }
 789
 790   if(strstr(option,"tim"))
 791     gBenchmark->Start("EMCALDigitizer");
 792
 793   AliRunLoader *rl = AliRunLoader::Instance();
 794
 795   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 796   if(!emcalLoader)
 797   {
 798     AliFatal("Did not get the  Loader");
 799     return; // coverity
 800   }
 801
 802   if (fLastEvent == -1)
 803     fLastEvent = rl->GetNumberOfEvents() - 1 ;
 804   else if (fDigInput)
 805     fLastEvent = fFirstEvent ; // what is this ??
 806
 807   Int_t nEvents = fLastEvent - fFirstEvent + 1;
 808   Int_t ievent  = -1;
 809
 810   AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
 811   if(!emcal)
 812   {
 813     AliFatal("Did not get the AliEMCAL pointer");
 814     return; // coverity
 815   }
 816
 817   AliEMCALGeometry *geom = emcal->GetGeometry();
 818   if(!geom)
 819   {
 820     AliFatal("Geometry pointer null");
 821     return; // fix for coverity
 822   }
 823
 824   const Int_t nTRU = geom->GetNTotalTRU();
 825   TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit",    nTRU*96);
 826   TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
 827
 828   rl->LoadSDigits("EMCAL");
 829   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++)
 830   {
 831     rl->GetEvent(ievent);
 832
 833     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
 834
 835     WriteDigits() ;
 836
 837     //Trigger Digits
 838     //-------------------------------------
 839
 840     Digits2FastOR(digitsTMP, digitsTRG);
 841
 842     WriteDigits(digitsTRG);
 843
 844     (emcalLoader->TreeD())->Fill();
 845
 846     emcalLoader->WriteDigits(   "OVERWRITE");
 847
 848     Unload();
 849
 850     digitsTRG  ->Delete();
 851     digitsTMP  ->Delete();
 852
 853     //-------------------------------------
 854
 855     if(strstr(option,"deb"))
 856       PrintDigits(option);
 857     if(strstr(option,"table")) gObjectTable->Print();
 858
 859     //increment the total number of Digits per run
 860     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
 861   }//loop
 862
 863   if(strstr(option,"tim"))
 864   {
 865     gBenchmark->Stop("EMCALDigitizer");
 866     Float_t cputime   = gBenchmark->GetCpuTime("EMCALDigitizer");
 867     Float_t avcputime = cputime;
 868     if(nEvents==0) avcputime = 0 ;
 869     AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
 870   }
 871 }
 872
 873 //__________________________________________________________________
 874 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
 875 {
 876   // Assign a smeared time to the digit, from observed distribution in LED system (?).
 877   // From F. Blanco
 878   Float_t res = -1;
 879   if (energy > 0)
 880   {
 881     res = TMath::Sqrt(fTimeResolutionPar0 +
 882                       fTimeResolutionPar1/(energy*energy) );
 883   }
 884   // parametrization above is for ns. Convert to seconds:
 885   return res*1e-9;
 886 }
 887
 888 //____________________________________________________________________________
 889 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
 890 {
 891   // FEE digits afterburner to produce TRG digits
 892   // we are only interested in the FEE digit deposited energy
 893   // to be converted later into a voltage value
 894
 895   // push the FEE digit to its associated FastOR (numbered from 0:95)
 896   // TRU is in charge of summing module digits
 897
 898   AliRunLoader *runLoader = AliRunLoader::Instance();
 899
 900   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
 901   if (!emcalLoader) AliFatal("Did not get the  Loader");
 902
 903   const  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 904   if(!geom)
 905   {
 906     AliFatal("Geometry pointer null");
 907     return; // fix for coverity
 908   }
 909
 910   // build FOR from simulated digits
 911   // and xfer to the corresponding TRU input (mapping)
 912
 913   TClonesArray* sdigits = emcalLoader->SDigits();
 914
 915   TClonesArray *digits = (TClonesArray*)sdigits->Clone();
 916
 917   AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
 918
 919   TIter NextDigit(digits);
 920   while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
 921   {
 922     if (IsDead(digit)) continue;
 923
 924     DecalibrateTrigger(digit);
 925
 926     Int_t id = digit->GetId();
 927
 928     Int_t trgid;
 929     if (geom->GetFastORIndexFromCellIndex(id , trgid))
 930     {
 931       AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
 932
 933       AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
 934
 935       if (!d)
 936       {
 937         new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
 938         d = (AliEMCALDigit*)digitsTMP->At(trgid);
 939         d->SetId(trgid);
 940       }
 941       else
 942       {
 943         *d = *d + *digit;
 944       }
 945     }
 946   }
 947
 948   if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
 949
 950   Int_t     nSamples = geom->GetNTotalTRU();
 951   Int_t *timeSamples = new Int_t[nSamples];
 952
 953   NextDigit = TIter(digitsTMP);
 954   while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
 955   {
 956     if (digit)
 957     {
 958       Int_t     id = digit->GetId();
 959       Float_t time = 50.e-9;
 960
 961       Double_t depositedEnergy = 0.;
 962       for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
 963
 964       if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
 965
 966       // FIXME: Check digit time!
 967       if (depositedEnergy) {
 968         depositedEnergy += gRandom->Gaus(0., .08);
 969         DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
 970
 971         for (Int_t j=0;j<nSamples;j++) {
 972           if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
 973           timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
 974         }
 975
 976         new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
 977
 978         if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
 979       }
 980     }
 981   }
 982
 983   delete [] timeSamples;
 984
 985   if (digits && digits->GetEntriesFast()) digits->Delete();
 986 }
 987
 988 //____________________________________________________________________________
 989 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
 990 {
 991   // parameters:
 992   // id: 0..95
 993   const Int_t    reso = 12;      // 11-bit resolution ADC
 994   const Double_t vFSR = 2.;      // Full scale input voltage range 2V (p-p)
 995 //const Double_t dNe  = 125;     // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
 996   const Double_t dNe  = 125/1.3; // F-ALTRO max V. FEE: factor 4
 997   const Double_t vA   = .136e-6; // CSP output range: 0.136uV/e-
 998   const Double_t rise = 50e-9;   // rise time (10-90%) of the FastOR signal before shaping
 999
1000   const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
1001
1002   Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
1003
1004   TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
1005   signalF.SetParameter( 0,   vV );
1006   signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
1007   signalF.SetParameter( 2, rise );
1008
1009   for (Int_t iTime=0; iTime<nSamples; iTime++)
1010   {
1011     // FIXME: add noise (probably not simply Gaussian) according to DA measurements
1012     // probably plan an access to OCDB
1013     Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
1014     if (TMath::Abs(sig) > vFSR/2.) {
1015       AliError("Signal overflow!");
1016       timeSamples[iTime] = (1 << reso) - 1;
1017     } else {
1018       AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
1019       timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
1020     }
1021   }
1022 }
1023
1024 //____________________________________________________________________________
1025 Bool_t AliEMCALDigitizer::Init()
1026 {
1027   // Makes all memory allocations
1028   fInit = kTRUE ;
1029   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1030
1031   if ( emcalLoader == 0 ) {
1032     Fatal("Init", "Could not obtain the AliEMCALLoader");
1033     return kFALSE;
1034   }
1035
1036   fFirstEvent = 0 ;
1037   fLastEvent = fFirstEvent ;
1038
1039   if(fDigInput)
1040     fInput = fDigInput->GetNinputs() ;
1041   else
1042     fInput           = 1 ;
1043
1044   fInputFileNames  = new TString[fInput] ;
1045   fEventNames      = new TString[fInput] ;
1046   fInputFileNames[0] = GetTitle() ;
1047   fEventNames[0]     = fEventFolderName.Data() ;
1048   Int_t index ;
1049   for (index = 1 ; index < fInput ; index++) {
1050     fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1051     TString tempo = fDigInput->GetInputFolderName(index) ;
1052     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
1053   }
1054
1055   //Calibration instance
1056   fCalibData = emcalLoader->CalibData();
1057   return fInit ;
1058 }
1059
1060 //____________________________________________________________________________
1061 void AliEMCALDigitizer::InitParameters()
1062 {
1063   // Parameter initialization for digitizer
1064
1065   // Get the parameters from the OCDB via the loader
1066   AliRunLoader *rl = AliRunLoader::Instance();
1067   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1068   AliEMCALSimParam * simParam = 0x0;
1069   if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1070
1071   if(!simParam){
1072     simParam = AliEMCALSimParam::GetInstance();
1073     AliWarning("Simulation Parameters not available in OCDB?");
1074   }
1075
1076   fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400;  // electrons per GeV
1077   fGainFluctuations   = simParam->GetGainFluctuations()   ; // 15.0;
1078
1079   fPinNoise           = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1080   if (fPinNoise < 0.0001 )
1081     Warning("InitParameters", "No noise added\n") ;
1082   fTimeNoise          = simParam->GetTimeNoise(); // 1.28E-5 s
1083   fDigitThreshold     = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1084   fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1085   fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1086   fTimeDelay          = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1087
1088   // These defaults are normally not used.
1089   // Values are read from calibration database instead
1090   fADCchannelEC       = 0.0162; // Nominal value set online for most of the channels, MIP peak sitting at ~266./16/1024
1091   fADCpedestalEC      = 0.0 ;   // GeV
1092   fADCchannelECDecal  = 1.0;    // No decalibration by default
1093   fTimeChannel        = 0.0;    // No time calibration by default
1094   fTimeChannelDecal   = 0.0;    // No time decalibration by default
1095
1096   fNADCEC             = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
1097
1098   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",
1099                   fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1100
1101 }
1102
1103 //__________________________________________________________________
1104 void AliEMCALDigitizer::Print1(Option_t * option)
1105 { // 19-nov-04 - just for convenience
1106   Print();
1107   PrintDigits(option);
1108 }
1109
1110 //__________________________________________________________________
1111 void AliEMCALDigitizer::Print (Option_t * ) const
1112 {
1113   // Print Digitizer's parameters
1114   printf("Print: \n------------------- %s -------------", GetName() ) ;
1115   if( strcmp(fEventFolderName.Data(), "") != 0 ){
1116     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
1117
1118     Int_t nStreams ;
1119     if (fDigInput)
1120       nStreams =  GetNInputStreams() ;
1121     else
1122       nStreams = fInput ;
1123
1124     AliRunLoader *rl=0;
1125
1126     Int_t index = 0 ;
1127     for (index = 0 ; index < nStreams ; index++) {
1128       TString tempo(fEventNames[index]) ;
1129       tempo += index ;
1130       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1131         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1132       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1133       if(emcalLoader){
1134         TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1135         if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1136           fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
1137         printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1138       }//loader
1139     }
1140
1141     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1142
1143     if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1144     else printf("\nNULL LOADER");
1145
1146     printf("\nWith following parameters:\n") ;
1147     printf("    Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1148     printf("    Threshold  in Tower  (fDigitThreshold) = %d\n", fDigitThreshold)  ;
1149     printf("---------------------------------------------------\n")  ;
1150   }
1151   else
1152     printf("Print: AliEMCALDigitizer not initialized") ;
1153 }
1154
1155 //__________________________________________________________________
1156 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1157 {
1158   //utility method for printing digit information
1159
1160   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1161   if(emcalLoader){
1162     TClonesArray * digits  = emcalLoader->Digits() ;
1163     TClonesArray * sdigits = emcalLoader->SDigits() ;
1164
1165     printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1166     printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1167
1168     if(strstr(option,"all")){
1169       //loop over digits
1170       AliEMCALDigit * digit;
1171       printf("\nEMC digits (with primaries):\n")  ;
1172       printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
1173       Int_t index ;
1174       for (index = 0 ; index < digits->GetEntries() ; index++) {
1175         digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1176         if(digit){
1177           printf("\n%6d  %8f    %6.5e %4d      %2d : ",
1178                  digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1179           Int_t iprimary;
1180           for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1181             printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1182           }
1183         }// digit exists
1184       }// loop
1185     }
1186     printf("\n");
1187   }// loader exists
1188   else printf("NULL LOADER, cannot print\n");
1189 }
1190
1191 //__________________________________________________________________
1192 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1193 {
1194   // Calculates the time signal generated by noise
1195   //printf("Time noise %e\n",fTimeNoise);
1196   return gRandom->Rndm() * fTimeNoise;
1197 }
1198
1199 //__________________________________________________________________
1200 void AliEMCALDigitizer::Unload()
1201 {
1202   // Unloads the SDigits and Digits
1203   AliRunLoader *rl=0;
1204
1205   Int_t i ;
1206   for(i = 1 ; i < fInput ; i++){
1207     TString tempo(fEventNames[i]) ;
1208     tempo += i;
1209     if ((rl = AliRunLoader::GetRunLoader(tempo)))
1210       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1211   }
1212   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1213   if(emcalLoader)emcalLoader->UnloadDigits() ;
1214   else AliFatal("Did not get the loader");
1215 }
1216
1217 //_________________________________________________________________________________________
1218 void AliEMCALDigitizer::WriteDigits()
1219 { // Makes TreeD in the output file.
1220   // Check if branch already exists:
1221   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1222   //      already existing branches.
1223   //   else creates branch with Digits, named "EMCAL", title "...",
1224   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1225   //      and names of files, from which digits are made.
1226
1227   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1228
1229   if(emcalLoader){
1230     const TClonesArray * digits = emcalLoader->Digits() ;
1231     TTree * treeD = emcalLoader->TreeD();
1232     if ( !treeD ) {
1233       emcalLoader->MakeDigitsContainer();
1234       treeD = emcalLoader->TreeD();
1235     }
1236
1237     // -- create Digits branch
1238     Int_t bufferSize = 32000 ;
1239     TBranch * digitsBranch = 0;
1240     if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1241       digitsBranch->SetAddress(&digits);
1242       AliWarning("Digits Branch already exists. Not all digits will be visible");
1243     }
1244     else
1245       treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1246     //digitsBranch->SetTitle(fEventFolderName);
1247
1248     //  treeD->Fill() ;
1249     /*
1250      emcalLoader->WriteDigits("OVERWRITE");
1251      emcalLoader->WriteDigitizer("OVERWRITE");
1252
1253      Unload() ;
1254      */
1255
1256   }// loader exists
1257   else AliFatal("Loader not available");
1258 }
1259
1260 //__________________________________________________________________
1261 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1262 {
1263   // overloaded method
1264   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1265   if(emcalLoader){
1266
1267     TTree* treeD = emcalLoader->TreeD();
1268     if (!treeD)
1269       {
1270         emcalLoader->MakeDigitsContainer();
1271         treeD = emcalLoader->TreeD();
1272       }
1273
1274     // -- create Digits branch
1275     Int_t bufferSize = 32000;
1276
1277     if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1278       {
1279         triggerBranch->SetAddress(&digits);
1280       }
1281     else
1282       {
1283         treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1284       }
1285
1286     //  treeD->Fill();
1287   }// loader exists
1288   else AliFatal("Loader not available");
1289 }
1290
1291 //__________________________________________________________________
1292 Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1293 {
1294   // Check if cell is defined as dead, so that it is not included
1295   // input is digit
1296   Int_t absId   = digit->GetId();
1297
1298   return IsDead(absId);
1299
1300 }
1301
1302
1303 //__________________________________________________________________
1304 Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1305 {
1306   // Check if cell absID is defined as dead, so that it is not included
1307
1308     AliRunLoader *runLoader = AliRunLoader::Instance();
1309     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1310     if (!emcalLoader)
1311     {
1312       AliFatal("Did not get the  Loader");
1313       return kTRUE;
1314     }
1315
1316     AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1317     if (!caloPed)
1318     {
1319         AliWarning("Could not access pedestal data! No dead channel removal applied");
1320         return kFALSE;
1321     }
1322
1323         // Load Geometry
1324     const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1325     if (!geom)
1326     {
1327       AliFatal("Did not get geometry from EMCALLoader");
1328       return kTRUE; //coverity
1329     }
1330
1331     Int_t iSupMod = -1;
1332     Int_t nModule = -1;
1333     Int_t nIphi   = -1;
1334     Int_t nIeta   = -1;
1335     Int_t iphi    = -1;
1336     Int_t ieta    = -1;
1337
1338     Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1339
1340     if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1341     geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1342
1343     Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1344
1345     if (channelStatus == AliCaloCalibPedestal::kDead)
1346         return kTRUE;
1347     else
1348         return kFALSE;
1349 }
1350