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