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