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