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