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