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