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