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