Adding include path to allow compilation of CleanGeom task
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigitizer.cxx
CommitLineData
61e0abb5 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//_________________________________________________________________________
ffa6d63b 19//
61e0abb5 20//////////////////////////////////////////////////////////////////////////////
ab6a174f 21// Class performs digitization of Summable digits from simulated data
ffa6d63b 22//
6cc75819 23// In addition it performs mixing/embedding of summable digits from different events.
61e0abb5 24//
61d80ce0 25// For each event 3 branches are created in TreeD:
61e0abb5 26// "EMCAL" - list of digits
61d80ce0 27// "EMCALTRG" - list of trigger digits
61e0abb5 28// "AliEMCALDigitizer" - AliEMCALDigitizer with all parameters used in digitization
29//
61e0abb5 30//
ffa6d63b 31////////////////////////////////////////////////////////////////////////////////////
61e0abb5 32//
ffa6d63b 33//*-- Author: Sahal Yacoob (LBL)
814ad4bf 34// based on : AliEMCALDigitizer
05a92d59 35// Modif:
36// August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
6cc75819 37// of new IO (a la PHOS)
1963b290 38// November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry
61d80ce0 39// July 2011 GCB: Digitizer modified to accomodate embedding.
40// Time calibration added. Decalibration possibility of energy and time added
1963b290 41//_________________________________________________________________________________
61e0abb5 42
43// --- ROOT system ---
e939a978 44#include <TROOT.h>
45#include <TTree.h>
46#include <TSystem.h>
47#include <TBenchmark.h>
e939a978 48#include <TBrowser.h>
49#include <TObjectTable.h>
50#include <TRandom.h>
916f1e76 51#include <TF1.h>
21e7df44 52#include <cassert>
61e0abb5 53
54// --- AliRoot header files ---
b42d4197 55#include "AliLog.h"
e4dcc484 56#include "AliRun.h"
f21fc003 57#include "AliDigitizationInput.h"
5dee926e 58#include "AliRunLoader.h"
53f1c378 59#include "AliCDBManager.h"
60#include "AliCDBEntry.h"
61e0abb5 61#include "AliEMCALDigit.h"
05a92d59 62#include "AliEMCAL.h"
5dee926e 63#include "AliEMCALLoader.h"
61e0abb5 64#include "AliEMCALDigitizer.h"
65#include "AliEMCALSDigitizer.h"
66#include "AliEMCALGeometry.h"
a435f763 67#include "AliEMCALCalibData.h"
33a52e55 68#include "AliEMCALSimParam.h"
916f1e76 69#include "AliEMCALRawDigit.h"
6995e8b7 70#include "AliCaloCalibPedestal.h"
916f1e76 71
61d80ce0 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 }
05a92d59 108
61e0abb5 109ClassImp(AliEMCALDigitizer)
61d80ce0 110
111
61e0abb5 112//____________________________________________________________________________
18a21c7c 113AliEMCALDigitizer::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),
5fb45456 123 fGainFluctuations(0),
18a21c7c 124 fPinNoise(0),
6cc75819 125 fTimeNoise(0),
f0a6dc6f 126 fTimeDelay(0),
a2f8e711 127 fTimeResolutionPar0(0),
128 fTimeResolutionPar1(0),
18a21c7c 129 fADCchannelEC(0),
130 fADCpedestalEC(0),
6cc75819 131 fADCchannelECDecal(0),
132 fTimeChannel(0),
133 fTimeChannelDecal(0),
18a21c7c 134 fNADCEC(0),
135 fEventFolderName(""),
136 fFirstEvent(0),
137 fLastEvent(0),
f21fc003 138 fCalibData(0x0),
139 fSDigitizer(0x0)
61e0abb5 140{
141 // ctor
88cb7938 142 InitParameters() ;
f21fc003 143 fDigInput = 0 ; // We work in the standalone mode
839828a6 144}
61e0abb5 145
839828a6 146//____________________________________________________________________________
18a21c7c 147AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
f21fc003 148 : AliDigitizer("EMCALDigitizer", alirunFileName),
18a21c7c 149 fDefaultInit(kFALSE),
150 fDigitsInRun(0),
151 fInit(0),
152 fInput(0),
153 fInputFileNames(0),
154 fEventNames(0),
155 fDigitThreshold(0),
156 fMeanPhotonElectron(0),
5fb45456 157 fGainFluctuations(0),
18a21c7c 158 fPinNoise(0),
6cc75819 159 fTimeNoise(0),
96d5ea0d 160 fTimeDelay(0),
a2f8e711 161 fTimeResolutionPar0(0),
162 fTimeResolutionPar1(0),
18a21c7c 163 fADCchannelEC(0),
164 fADCpedestalEC(0),
6cc75819 165 fADCchannelECDecal(0),
166 fTimeChannel(0),
167 fTimeChannelDecal(0),
18a21c7c 168 fNADCEC(0),
169 fEventFolderName(eventFolderName),
170 fFirstEvent(0),
171 fLastEvent(0),
f21fc003 172 fCalibData(0x0),
173 fSDigitizer(0x0)
839828a6 174{
05a92d59 175 // ctor
a9485c37 176 InitParameters() ;
88cb7938 177 Init() ;
f21fc003 178 fDigInput = 0 ; // We work in the standalone mode
61e0abb5 179}
839828a6 180
61e0abb5 181//____________________________________________________________________________
18a21c7c 182AliEMCALDigitizer::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),
5fb45456 192 fGainFluctuations(d.fGainFluctuations),
18a21c7c 193 fPinNoise(d.fPinNoise),
6cc75819 194 fTimeNoise(d.fTimeNoise),
f0a6dc6f 195 fTimeDelay(d.fTimeDelay),
a2f8e711 196 fTimeResolutionPar0(d.fTimeResolutionPar0),
197 fTimeResolutionPar1(d.fTimeResolutionPar1),
18a21c7c 198 fADCchannelEC(d.fADCchannelEC),
199 fADCpedestalEC(d.fADCpedestalEC),
6cc75819 200 fADCchannelECDecal(d.fADCchannelECDecal),
201 fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
18a21c7c 202 fNADCEC(d.fNADCEC),
203 fEventFolderName(d.fEventFolderName),
204 fFirstEvent(d.fFirstEvent),
205 fLastEvent(d.fLastEvent),
f21fc003 206 fCalibData(d.fCalibData),
207 fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
88cb7938 208{
209 // copyy ctor
96d5ea0d 210}
88cb7938 211
212//____________________________________________________________________________
f21fc003 213AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
214 : AliDigitizer(rd,"EMCALDigitizer"),
18a21c7c 215 fDefaultInit(kFALSE),
216 fDigitsInRun(0),
217 fInit(0),
218 fInput(0),
219 fInputFileNames(0),
220 fEventNames(0),
da509d54 221 fDigitThreshold(0),
18a21c7c 222 fMeanPhotonElectron(0),
5fb45456 223 fGainFluctuations(0),
da509d54 224 fPinNoise(0.),
6cc75819 225 fTimeNoise(0.),
f0a6dc6f 226 fTimeDelay(0.),
a2f8e711 227 fTimeResolutionPar0(0.),
228 fTimeResolutionPar1(0.),
18a21c7c 229 fADCchannelEC(0),
230 fADCpedestalEC(0),
6cc75819 231 fADCchannelECDecal(0),
232 fTimeChannel(0),
233 fTimeChannelDecal(0),
18a21c7c 234 fNADCEC(0),
235 fEventFolderName(0),
236 fFirstEvent(0),
237 fLastEvent(0),
f21fc003 238 fCalibData(0x0),
239 fSDigitizer(0x0)
61e0abb5 240{
45fa49ca 241 // ctor Init() is called by RunDigitizer
f21fc003 242 fDigInput = rd ;
243 fEventFolderName = fDigInput->GetInputFolderName(0) ;
244 SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
05a92d59 245 InitParameters() ;
839828a6 246}
814ad4bf 247
839828a6 248//____________________________________________________________________________
249 AliEMCALDigitizer::~AliEMCALDigitizer()
250{
14ce0a6e 251 //dtor
f5f384f4 252 delete [] fInputFileNames ;
253 delete [] fEventNames ;
f21fc003 254 if (fSDigitizer) delete fSDigitizer;
61e0abb5 255}
256
257//____________________________________________________________________________
a81281f7 258void AliEMCALDigitizer::Digitize(Int_t event)
259{
61e0abb5 260 // Makes the digitization of the collected summable digits
261 // for this it first creates the array of all EMCAL modules
a81281f7 262 // filled with noise and after that adds contributions from
263 // SDigits. This design helps to avoid scanning over the
ab6a174f 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 //
f69e318a 271
33c3c91a 272 AliRunLoader *rl = AliRunLoader::Instance();
5dee926e 273 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
7e1d9a9b 274
f69e318a 275 if(!emcalLoader)
a81281f7 276 {
f69e318a 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
a81281f7 283 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
f69e318a 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
cd9ed873 310 //-------------------------------------------------------
f69e318a 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 ;
6cc75819 321
f69e318a 322 AliRunLoader * rl2 = AliRunLoader::GetRunLoader(tempo) ;
6cc75819 323
cd9ed873 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 }
96d5ea0d 332
f69e318a 333 if (fDigInput)
cd9ed873 334 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
335
f69e318a 336 Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
cd9ed873 337
f69e318a 338 rl2->LoadSDigits();
339 // rl2->LoadDigits();
340 rl2->GetEvent(readEvent);
cd9ed873 341
342 if(!rl2->GetDetectorLoader("EMCAL"))
a81281f7 343 {
cd9ed873 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 ...
f69e318a 354 }
cd9ed873 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
f69e318a 370 }// input loop
371
cd9ed873 372 //--------------------------------
f69e318a 373 //Find the first tower with signal
374 Int_t nextSig = nEMC + 1 ;
375 TClonesArray * sdigits ;
cd9ed873 376 for(i = 0 ; i < fInput ; i++)
377 {
f69e318a 378 if(i > 0 && embed) continue;
8e708def 379
f69e318a 380 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
8e708def 381 if(!sdigits)
f69e318a 382 {
383 AliDebug(1,"Null sdigit pointer");
384 continue;
385 }
8e708def 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
f69e318a 405 }// input loop
8e708def 406
f69e318a 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
cd9ed873 415 //---------------------------------------------
f69e318a 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
a81281f7 421
f69e318a 422 if (IsDead(absID)) continue; // Don't instantiate dead digits
a81281f7 423
f69e318a 424 Float_t energy = 0 ;
a81281f7 425
f69e318a 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
8e708def 433 if (!digit)
434 {
435 AliDebug(1,"Digit pointer is null");
436 continue;
437 }
438
439 if(absID==nextSig)
f69e318a 440 {
8e708def 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++)
a81281f7 449 {
8e708def 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 ;
f69e318a 458
8e708def 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)
f69e318a 468 {
8e708def 469 aTime = curSDigit->GetAmplitude();
470 time = curSDigit->GetTime();
f69e318a 471 }
8e708def 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
f69e318a 482
8e708def 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 ;
f69e318a 490
8e708def 491 //calculate and set time
492 digit->SetTime(time) ;
6cc75819 493
8e708def 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]));
83366f85 506 if ( tmpdigit ) curNext = tmpdigit->GetId() ;
8e708def 507 }
83366f85 508
8e708def 509 if(curNext < nextSig) nextSig = curNext ;
510 }// sdigits exist
511 } // input loop
f69e318a 512
8e708def 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) ;
a81281f7 532
f69e318a 533 } // for(absID = 0; absID < nEMC; absID++)
534
cd9ed873 535 //---------------------------------------------------------
f69e318a 536 //Embed simulated amplitude (and time?) to real data digits
83366f85 537 if ( embed )
f69e318a 538 {
539 for(Int_t input = 1; input<fInput; input++)
a81281f7 540 {
f69e318a 541 TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
542 if(!realDigits)
a81281f7 543 {
f69e318a 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) ) ;
83366f85 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()));
f69e318a 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
cd9ed873 589 //------------------------------
f69e318a 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) ) ;
83366f85 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)
f69e318a 609 digits->RemoveAt(i) ;
83366f85 610
f69e318a 611 } // digit loop
612
613 digits->Compress() ;
614
615 Int_t ndigits = digits->GetEntriesFast() ;
616
cd9ed873 617 //---------------------------------------------------------------
f69e318a 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) ) ;
83366f85 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());
f69e318a 639 }//Digit loop
6cc75819 640
814ad4bf 641}
61e0abb5 642
6cc75819 643//_____________________________________________________________________
644void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
a81281f7 645{
ab6a174f 646 // JLK 26-June-2008
e4dcc484 647 // Returns digitized value of the energy in a cell absId
ab6a174f 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 //
53f1c378 653 // Load Geometry
654 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
d62a891a 655
a81281f7 656 if (geom==0)
657 {
96d5ea0d 658 AliFatal("Did not get geometry from EMCALLoader");
6cc75819 659 return;
e4dcc484 660 }
5eb74a51 661
6cc75819 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)
a81281f7 670 Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
6cc75819 671 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
96d5ea0d 672
a81281f7 673 if(fCalibData)
674 {
6cc75819 675 fADCpedestalEC = fCalibData->GetADCpedestal (iSupMod,ieta,iphi);
676 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
677 fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
031c3928 678 fTimeChannel = fCalibData->GetTimeChannel (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
a81281f7 679 fTimeChannelDecal = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
6cc75819 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
6337857a 686 if ( energy > fNADCEC ) energy = fNADCEC ;
6cc75819 687}
829ba234 688
6cc75819 689//_____________________________________________________________________
6337857a 690void AliEMCALDigitizer::DecalibrateTrigger(AliEMCALDigit *digit)
a81281f7 691{
692 // Decalibrate, used in Trigger digits
6337857a 693
694 if ( !fCalibData ) return ;
695
a81281f7 696 // Load Geometry
697 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
63c22917 698
a81281f7 699 if (!geom)
700 {
701 AliFatal("Did not get geometry from EMCALLoader");
702 return;
703 }
63c22917 704
6337857a 705 // Recover the digit location in row-column-SM
a81281f7 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;
52e22b66 713
a81281f7 714 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
a81281f7 715 if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
6337857a 716
a81281f7 717 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
63c22917 718
6337857a 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
63c22917 727}
728
729//_____________________________________________________________________
6cc75819 730void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
a81281f7 731{
6cc75819 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
a81281f7 743 if (!geom)
744 {
6cc75819 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) ;
a81281f7 756 if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
6cc75819 757 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
758
a81281f7 759 if(fCalibData)
760 {
6cc75819 761 fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
762 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
031c3928 763 fTimeChannel = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
6cc75819 764 }
765
766 adc = adc * fADCchannelEC - fADCpedestalEC;
767 time -= fTimeChannel;
96d5ea0d 768
61e0abb5 769}
770
6cc75819 771
61e0abb5 772//____________________________________________________________________________
2f9aea1a 773void AliEMCALDigitizer::Digitize(Option_t *option)
774{
45fa49ca 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)
2f9aea1a 780
781 if (!fInit)
782 { // to prevent overwrite existing file
f21fc003 783 Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
88cb7938 784 return ;
2f9aea1a 785 }
786
787 if (strstr(option,"print"))
788 {
88cb7938 789 Print();
2f9aea1a 790 return ;
814ad4bf 791 }
792
61e0abb5 793 if(strstr(option,"tim"))
794 gBenchmark->Start("EMCALDigitizer");
2f9aea1a 795
33c3c91a 796 AliRunLoader *rl = AliRunLoader::Instance();
2f9aea1a 797
5dee926e 798 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
efb8df82 799 if(!emcalLoader)
800 {
96d5ea0d 801 AliFatal("Did not get the Loader");
2f9aea1a 802 return; // coverity
30bc5594 803 }
2f9aea1a 804
805 if (fLastEvent == -1)
806 fLastEvent = rl->GetNumberOfEvents() - 1 ;
807 else if (fDigInput)
808 fLastEvent = fFirstEvent ; // what is this ??
809
810 Int_t nEvents = fLastEvent - fFirstEvent + 1;
811 Int_t ievent = -1;
812
813 AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
814 if(!emcal)
efb8df82 815 {
2f9aea1a 816 AliFatal("Did not get the AliEMCAL pointer");
817 return; // coverity
818 }
819
820 AliEMCALGeometry *geom = emcal->GetGeometry();
821 if(!geom)
822 {
823 AliFatal("Geometry pointer null");
824 return; // fix for coverity
825 }
826
827 const Int_t nTRU = geom->GetNTotalTRU();
828 TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit", nTRU*96);
829 TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
830
831 rl->LoadSDigits("EMCAL");
832 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++)
833 {
834 rl->GetEvent(ievent);
96d5ea0d 835
2f9aea1a 836 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
efb8df82 837
2f9aea1a 838 WriteDigits() ;
839
840 //Trigger Digits
841 //-------------------------------------
842
843 Digits2FastOR(digitsTMP, digitsTRG);
844
845 WriteDigits(digitsTRG);
846
847 (emcalLoader->TreeD())->Fill();
848
849 emcalLoader->WriteDigits( "OVERWRITE");
850
851 Unload();
852
853 digitsTRG ->Delete();
854 digitsTMP ->Delete();
855
856 //-------------------------------------
857
858 if(strstr(option,"deb"))
859 PrintDigits(option);
860 if(strstr(option,"table")) gObjectTable->Print();
861
862 //increment the total number of Digits per run
863 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
864 }//loop
64942713 865
2f9aea1a 866 if(strstr(option,"tim"))
867 {
61e0abb5 868 gBenchmark->Stop("EMCALDigitizer");
0bd264d5 869 Float_t cputime = gBenchmark->GetCpuTime("EMCALDigitizer");
870 Float_t avcputime = cputime;
871 if(nEvents==0) avcputime = 0 ;
872 AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
2f9aea1a 873 }
61e0abb5 874}
875
a2f8e711 876//__________________________________________________________________
877Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
878{
a81281f7 879 // Assign a smeared time to the digit, from observed distribution in LED system (?).
a2f8e711 880 // From F. Blanco
881 Float_t res = -1;
a81281f7 882 if (energy > 0)
883 {
a2f8e711 884 res = TMath::Sqrt(fTimeResolutionPar0 +
885 fTimeResolutionPar1/(energy*energy) );
886 }
b2a891e0 887 // parametrization above is for ns. Convert to seconds:
888 return res*1e-9;
a2f8e711 889}
890
a81281f7 891//____________________________________________________________________________
916f1e76 892void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
893{
a81281f7 894 // FEE digits afterburner to produce TRG digits
de39a0ff 895 // we are only interested in the FEE digit deposited energy
896 // to be converted later into a voltage value
897
898 // push the FEE digit to its associated FastOR (numbered from 0:95)
899 // TRU is in charge of summing module digits
900
901 AliRunLoader *runLoader = AliRunLoader::Instance();
902
de39a0ff 903 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
6995e8b7 904 if (!emcalLoader) AliFatal("Did not get the Loader");
905
a81281f7 906 const AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
efb8df82 907 if(!geom)
908 {
909 AliFatal("Geometry pointer null");
910 return; // fix for coverity
911 }
a81281f7 912
913 // build FOR from simulated digits
914 // and xfer to the corresponding TRU input (mapping)
915
916 TClonesArray* sdigits = emcalLoader->SDigits();
917
918 TClonesArray *digits = (TClonesArray*)sdigits->Clone();
919
920 AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
921
922 TIter NextDigit(digits);
923 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
924 {
925 if (IsDead(digit)) continue;
926
6337857a 927 DecalibrateTrigger(digit);
a81281f7 928
929 Int_t id = digit->GetId();
930
931 Int_t trgid;
efb8df82 932 if (geom->GetFastORIndexFromCellIndex(id , trgid))
a81281f7 933 {
934 AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
7e1d9a9b 935
a81281f7 936 AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
7d525d9d 937
a81281f7 938 if (!d)
96d5ea0d 939 {
a81281f7 940 new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
941 d = (AliEMCALDigit*)digitsTMP->At(trgid);
942 d->SetId(trgid);
943 }
944 else
945 {
946 *d = *d + *digit;
96d5ea0d 947 }
a81281f7 948 }
949 }
950
951 if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
952
953 Int_t nSamples = geom->GetNTotalTRU();
954 Int_t *timeSamples = new Int_t[nSamples];
955
956 NextDigit = TIter(digitsTMP);
957 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
958 {
959 if (digit)
960 {
961 Int_t id = digit->GetId();
962 Float_t time = 50.e-9;
7e1d9a9b 963
a81281f7 964 Double_t depositedEnergy = 0.;
965 for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
7e1d9a9b 966
a81281f7 967 if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
968
969 // FIXME: Check digit time!
970 if (depositedEnergy) {
971 depositedEnergy += gRandom->Gaus(0., .08);
972 DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
973
974 for (Int_t j=0;j<nSamples;j++) {
975 if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
976 timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
7e1d9a9b 977 }
a81281f7 978
979 new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
980
981 if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
96d5ea0d 982 }
a81281f7 983 }
984 }
985
986 delete [] timeSamples;
987
988 if (digits && digits->GetEntriesFast()) digits->Delete();
916f1e76 989}
990
a81281f7 991//____________________________________________________________________________
916f1e76 992void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
993{
a81281f7 994 // parameters:
995 // id: 0..95
996 const Int_t reso = 12; // 11-bit resolution ADC
997 const Double_t vFSR = 2.; // Full scale input voltage range 2V (p-p)
998//const Double_t dNe = 125; // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
999 const Double_t dNe = 125/1.3; // F-ALTRO max V. FEE: factor 4
1000 const Double_t vA = .136e-6; // CSP output range: 0.136uV/e-
1001 const Double_t rise = 50e-9; // rise time (10-90%) of the FastOR signal before shaping
916f1e76 1002
a81281f7 1003 const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
916f1e76 1004
a81281f7 1005 Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
1006
1007 TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
1008 signalF.SetParameter( 0, vV );
1009 signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
1010 signalF.SetParameter( 2, rise );
916f1e76 1011
a81281f7 1012 for (Int_t iTime=0; iTime<nSamples; iTime++)
1013 {
1014 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
1015 // probably plan an access to OCDB
28bcaaaa 1016 Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
1017 if (TMath::Abs(sig) > vFSR/2.) {
1018 AliError("Signal overflow!");
1019 timeSamples[iTime] = (1 << reso) - 1;
1020 } else {
1021 AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
1022 timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
1023 }
a81281f7 1024 }
916f1e76 1025}
1026
05a92d59 1027//____________________________________________________________________________
1028Bool_t AliEMCALDigitizer::Init()
1029{
1030 // Makes all memory allocations
88cb7938 1031 fInit = kTRUE ;
33c3c91a 1032 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
61d80ce0 1033
5dee926e 1034 if ( emcalLoader == 0 ) {
1035 Fatal("Init", "Could not obtain the AliEMCALLoader");
05a92d59 1036 return kFALSE;
1037 }
61d80ce0 1038
45fa49ca 1039 fFirstEvent = 0 ;
1040 fLastEvent = fFirstEvent ;
1041
f21fc003 1042 if(fDigInput)
1043 fInput = fDigInput->GetNinputs() ;
88cb7938 1044 else
1045 fInput = 1 ;
1046
1047 fInputFileNames = new TString[fInput] ;
1048 fEventNames = new TString[fInput] ;
1049 fInputFileNames[0] = GetTitle() ;
1050 fEventNames[0] = fEventFolderName.Data() ;
1051 Int_t index ;
1052 for (index = 1 ; index < fInput ; index++) {
f21fc003 1053 fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1054 TString tempo = fDigInput->GetInputFolderName(index) ;
1055 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
05a92d59 1056 }
f21fc003 1057
53f1c378 1058 //Calibration instance
1059 fCalibData = emcalLoader->CalibData();
88cb7938 1060 return fInit ;
05a92d59 1061}
1062
1063//____________________________________________________________________________
1064void AliEMCALDigitizer::InitParameters()
14ce0a6e 1065{
29b5d492 1066 // Parameter initialization for digitizer
6569f329 1067
1068 // Get the parameters from the OCDB via the loader
1069 AliRunLoader *rl = AliRunLoader::Instance();
1070 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1071 AliEMCALSimParam * simParam = 0x0;
1072 if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1073
1074 if(!simParam){
61d80ce0 1075 simParam = AliEMCALSimParam::GetInstance();
1076 AliWarning("Simulation Parameters not available in OCDB?");
6569f329 1077 }
61d80ce0 1078
5fb45456 1079 fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400; // electrons per GeV
1080 fGainFluctuations = simParam->GetGainFluctuations() ; // 15.0;
1081
6569f329 1082 fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
33a52e55 1083 if (fPinNoise < 0.0001 )
88cb7938 1084 Warning("InitParameters", "No noise added\n") ;
6cc75819 1085 fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
6569f329 1086 fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
a2f8e711 1087 fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1088 fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
f0a6dc6f 1089 fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
61d80ce0 1090
97075cd8 1091 // These defaults are normally not used.
1092 // Values are read from calibration database instead
6337857a 1093 fADCchannelEC = 0.0162; // Nominal value set online for most of the channels, MIP peak sitting at ~266./16/1024
6cc75819 1094 fADCpedestalEC = 0.0 ; // GeV
1095 fADCchannelECDecal = 1.0; // No decalibration by default
1096 fTimeChannel = 0.0; // No time calibration by default
1097 fTimeChannelDecal = 0.0; // No time decalibration by default
33a52e55 1098
6cc75819 1099 fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
61d80ce0 1100
5fb45456 1101 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",
1102 fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1103
05a92d59 1104}
61e0abb5 1105
61e0abb5 1106//__________________________________________________________________
1963b290 1107void AliEMCALDigitizer::Print1(Option_t * option)
6cc75819 1108{ // 19-nov-04 - just for convenience
1963b290 1109 Print();
1110 PrintDigits(option);
1111}
1112
61e0abb5 1113//__________________________________________________________________
6cc75819 1114void AliEMCALDigitizer::Print (Option_t * ) const
88cb7938 1115{
1116 // Print Digitizer's parameters
fdebddeb 1117 printf("Print: \n------------------- %s -------------", GetName() ) ;
88cb7938 1118 if( strcmp(fEventFolderName.Data(), "") != 0 ){
1119 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1120
1121 Int_t nStreams ;
f21fc003 1122 if (fDigInput)
88cb7938 1123 nStreams = GetNInputStreams() ;
1124 else
1125 nStreams = fInput ;
61d80ce0 1126
53e430a3 1127 AliRunLoader *rl=0;
88cb7938 1128
1129 Int_t index = 0 ;
1130 for (index = 0 ; index < nStreams ; index++) {
1131 TString tempo(fEventNames[index]) ;
1132 tempo += index ;
5dee926e 1133 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
96d5ea0d 1134 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
5dee926e 1135 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
96d5ea0d 1136 if(emcalLoader){
1137 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1138 if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1139 fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
1140 printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1141 }//loader
88cb7938 1142 }
61d80ce0 1143
33c3c91a 1144 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
61d80ce0 1145
96d5ea0d 1146 if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1147 else printf("\nNULL LOADER");
61d80ce0 1148
88cb7938 1149 printf("\nWith following parameters:\n") ;
6cc75819 1150 printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
4354827d 1151 printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
88cb7938 1152 printf("---------------------------------------------------\n") ;
ea5f3389 1153 }
1154 else
fdebddeb 1155 printf("Print: AliEMCALDigitizer not initialized") ;
61e0abb5 1156}
814ad4bf 1157
61e0abb5 1158//__________________________________________________________________
14ce0a6e 1159void AliEMCALDigitizer::PrintDigits(Option_t * option)
1160{
1161 //utility method for printing digit information
61d80ce0 1162
33c3c91a 1163 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
96d5ea0d 1164 if(emcalLoader){
1165 TClonesArray * digits = emcalLoader->Digits() ;
1166 TClonesArray * sdigits = emcalLoader->SDigits() ;
1167
1168 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1169 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1170
1171 if(strstr(option,"all")){
1172 //loop over digits
1173 AliEMCALDigit * digit;
1174 printf("\nEMC digits (with primaries):\n") ;
1175 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1176 Int_t index ;
1177 for (index = 0 ; index < digits->GetEntries() ; index++) {
1178 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1179 if(digit){
1180 printf("\n%6d %8f %6.5e %4d %2d : ",
1181 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1182 Int_t iprimary;
1183 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1184 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1185 }
1186 }// digit exists
1187 }// loop
1188 }
1189 printf("\n");
1190 }// loader exists
1191 else printf("NULL LOADER, cannot print\n");
61e0abb5 1192}
05a92d59 1193
1194//__________________________________________________________________
1195Float_t AliEMCALDigitizer::TimeOfNoise(void)
04f0bda3 1196{
1197 // Calculates the time signal generated by noise
6cc75819 1198 //printf("Time noise %e\n",fTimeNoise);
1199 return gRandom->Rndm() * fTimeNoise;
05a92d59 1200}
1201
88cb7938 1202//__________________________________________________________________
1203void AliEMCALDigitizer::Unload()
1204{
fdebddeb 1205 // Unloads the SDigits and Digits
5dee926e 1206 AliRunLoader *rl=0;
1207
88cb7938 1208 Int_t i ;
1209 for(i = 1 ; i < fInput ; i++){
1210 TString tempo(fEventNames[i]) ;
5dee926e 1211 tempo += i;
1212 if ((rl = AliRunLoader::GetRunLoader(tempo)))
1213 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
88cb7938 1214 }
33c3c91a 1215 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
96d5ea0d 1216 if(emcalLoader)emcalLoader->UnloadDigits() ;
1217 else AliFatal("Did not get the loader");
88cb7938 1218}
1219
814ad4bf 1220//_________________________________________________________________________________________
9e5d2067 1221void AliEMCALDigitizer::WriteDigits()
a2f8e711 1222{ // Makes TreeD in the output file.
814ad4bf 1223 // Check if branch already exists:
1224 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1225 // already existing branches.
1226 // else creates branch with Digits, named "EMCAL", title "...",
1227 // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1228 // and names of files, from which digits are made.
96d5ea0d 1229
33c3c91a 1230 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
96d5ea0d 1231
1232 if(emcalLoader){
1233 const TClonesArray * digits = emcalLoader->Digits() ;
1234 TTree * treeD = emcalLoader->TreeD();
1235 if ( !treeD ) {
1236 emcalLoader->MakeDigitsContainer();
1237 treeD = emcalLoader->TreeD();
1238 }
1239
1240 // -- create Digits branch
1241 Int_t bufferSize = 32000 ;
1242 TBranch * digitsBranch = 0;
1243 if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1244 digitsBranch->SetAddress(&digits);
1245 AliWarning("Digits Branch already exists. Not all digits will be visible");
1246 }
1247 else
1248 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1249 //digitsBranch->SetTitle(fEventFolderName);
1250
1251 // treeD->Fill() ;
1252 /*
1253 emcalLoader->WriteDigits("OVERWRITE");
1254 emcalLoader->WriteDigitizer("OVERWRITE");
1255
1256 Unload() ;
1257 */
1258
1259 }// loader exists
1260 else AliFatal("Loader not available");
916f1e76 1261}
88cb7938 1262
916f1e76 1263//__________________________________________________________________
1264void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
efb8df82 1265{
1266 // overloaded method
61d80ce0 1267 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1268 if(emcalLoader){
96d5ea0d 1269
1270 TTree* treeD = emcalLoader->TreeD();
1271 if (!treeD)
61d80ce0 1272 {
1273 emcalLoader->MakeDigitsContainer();
1274 treeD = emcalLoader->TreeD();
1275 }
96d5ea0d 1276
1277 // -- create Digits branch
1278 Int_t bufferSize = 32000;
1279
1280 if (TBranch* triggerBranch = treeD->GetBranch(branchName))
61d80ce0 1281 {
1282 triggerBranch->SetAddress(&digits);
1283 }
96d5ea0d 1284 else
61d80ce0 1285 {
1286 treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1287 }
96d5ea0d 1288
1289 // treeD->Fill();
1290 }// loader exists
1291 else AliFatal("Loader not available");
814ad4bf 1292}
61e0abb5 1293
6995e8b7 1294//__________________________________________________________________
1295Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1296{
efb8df82 1297 // Check if cell is defined as dead, so that it is not included
1298 // input is digit
6995e8b7 1299 Int_t absId = digit->GetId();
6995e8b7 1300
efb8df82 1301 return IsDead(absId);
1302
6995e8b7 1303}
1304
1305
1306//__________________________________________________________________
1307Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1308{
efb8df82 1309 // Check if cell absID is defined as dead, so that it is not included
1310
6995e8b7 1311 AliRunLoader *runLoader = AliRunLoader::Instance();
1312 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
efb8df82 1313 if (!emcalLoader)
1314 {
1315 AliFatal("Did not get the Loader");
1316 return kTRUE;
1317 }
1318
6995e8b7 1319 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
efb8df82 1320 if (!caloPed)
1321 {
6995e8b7 1322 AliWarning("Could not access pedestal data! No dead channel removal applied");
1323 return kFALSE;
1324 }
1325
1326 // Load Geometry
1327 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
efb8df82 1328 if (!geom)
1329 {
1330 AliFatal("Did not get geometry from EMCALLoader");
1331 return kTRUE; //coverity
1332 }
1333
6995e8b7 1334 Int_t iSupMod = -1;
1335 Int_t nModule = -1;
1336 Int_t nIphi = -1;
1337 Int_t nIeta = -1;
1338 Int_t iphi = -1;
1339 Int_t ieta = -1;
1340
1341 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1342
1343 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1344 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1345
1346 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1347
1348 if (channelStatus == AliCaloCalibPedestal::kDead)
1349 return kTRUE;
1350 else
1351 return kFALSE;
1352}
1353