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