track number corrected
[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>
21e7df44 69#include <cassert>
61e0abb5 70
71// --- AliRoot header files ---
b42d4197 72#include "AliLog.h"
e4dcc484 73#include "AliRun.h"
ea5f3389 74#include "AliRunDigitizer.h"
5dee926e 75#include "AliRunLoader.h"
53f1c378 76#include "AliCDBManager.h"
77#include "AliCDBEntry.h"
61e0abb5 78#include "AliEMCALDigit.h"
05a92d59 79#include "AliEMCAL.h"
5dee926e 80#include "AliEMCALLoader.h"
61e0abb5 81#include "AliEMCALDigitizer.h"
82#include "AliEMCALSDigitizer.h"
83#include "AliEMCALGeometry.h"
05a92d59 84#include "AliEMCALTick.h"
85
61e0abb5 86ClassImp(AliEMCALDigitizer)
87
88
89//____________________________________________________________________________
18a21c7c 90AliEMCALDigitizer::AliEMCALDigitizer()
91 : AliDigitizer("",""),
92 fDefaultInit(kTRUE),
93 fDigitsInRun(0),
94 fInit(0),
95 fInput(0),
96 fInputFileNames(0x0),
97 fEventNames(0x0),
98 fDigitThreshold(0),
99 fMeanPhotonElectron(0),
100 fPedestal(0),
101 fSlope(0),
102 fPinNoise(0),
103 fTimeResolution(0),
104 fTimeThreshold(0),
105 fTimeSignalLength(0),
106 fADCchannelEC(0),
107 fADCpedestalEC(0),
108 fNADCEC(0),
109 fEventFolderName(""),
110 fFirstEvent(0),
111 fLastEvent(0),
7ea6391b 112 fCalibData(0x0)
61e0abb5 113{
114 // ctor
88cb7938 115 InitParameters() ;
ab6a174f 116 fManager = 0 ; // We work in the standalone mode
839828a6 117}
61e0abb5 118
839828a6 119//____________________________________________________________________________
18a21c7c 120AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
121 : AliDigitizer("EMCAL"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
122 fDefaultInit(kFALSE),
123 fDigitsInRun(0),
124 fInit(0),
125 fInput(0),
126 fInputFileNames(0),
127 fEventNames(0),
128 fDigitThreshold(0),
129 fMeanPhotonElectron(0),
130 fPedestal(0),
131 fSlope(0),
132 fPinNoise(0),
133 fTimeResolution(0),
134 fTimeThreshold(0),
135 fTimeSignalLength(0),
136 fADCchannelEC(0),
137 fADCpedestalEC(0),
138 fNADCEC(0),
139 fEventFolderName(eventFolderName),
140 fFirstEvent(0),
141 fLastEvent(0),
7ea6391b 142 fCalibData(0x0)
839828a6 143{
05a92d59 144 // ctor
a9485c37 145 InitParameters() ;
88cb7938 146 Init() ;
ab6a174f 147 fManager = 0 ; // We work in the standalone mode
61e0abb5 148}
839828a6 149
61e0abb5 150//____________________________________________________________________________
18a21c7c 151AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
152 : AliDigitizer(d.GetName(),d.GetTitle()),
153 fDefaultInit(d.fDefaultInit),
154 fDigitsInRun(d.fDigitsInRun),
155 fInit(d.fInit),
156 fInput(d.fInput),
157 fInputFileNames(d.fInputFileNames),
158 fEventNames(d.fEventNames),
159 fDigitThreshold(d.fDigitThreshold),
160 fMeanPhotonElectron(d.fMeanPhotonElectron),
161 fPedestal(d.fPedestal),
162 fSlope(d.fSlope),
163 fPinNoise(d.fPinNoise),
164 fTimeResolution(d.fTimeResolution),
165 fTimeThreshold(d.fTimeThreshold),
166 fTimeSignalLength(d.fTimeSignalLength),
167 fADCchannelEC(d.fADCchannelEC),
168 fADCpedestalEC(d.fADCpedestalEC),
169 fNADCEC(d.fNADCEC),
170 fEventFolderName(d.fEventFolderName),
171 fFirstEvent(d.fFirstEvent),
172 fLastEvent(d.fLastEvent),
7ea6391b 173 fCalibData(d.fCalibData)
88cb7938 174{
175 // copyy ctor
88cb7938 176 }
177
178//____________________________________________________________________________
18a21c7c 179AliEMCALDigitizer::AliEMCALDigitizer(AliRunDigitizer * rd)
180 : AliDigitizer(rd,"EMCAL"+AliConfig::Instance()->GetDigitizerTaskName()),
181 fDefaultInit(kFALSE),
182 fDigitsInRun(0),
183 fInit(0),
184 fInput(0),
185 fInputFileNames(0),
186 fEventNames(0),
187 fDigitThreshold(0.),
188 fMeanPhotonElectron(0),
189 fPedestal(0),
190 fSlope(0.),
191 fPinNoise(0),
192 fTimeResolution(0.),
193 fTimeThreshold(0),
194 fTimeSignalLength(0),
195 fADCchannelEC(0),
196 fADCpedestalEC(0),
197 fNADCEC(0),
198 fEventFolderName(0),
199 fFirstEvent(0),
200 fLastEvent(0),
7ea6391b 201 fCalibData(0x0)
61e0abb5 202{
45fa49ca 203 // ctor Init() is called by RunDigitizer
88cb7938 204 fManager = rd ;
e0e18b6d 205 fEventFolderName = fManager->GetInputFolderName(0) ;
88cb7938 206 SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
05a92d59 207 InitParameters() ;
839828a6 208}
814ad4bf 209
839828a6 210//____________________________________________________________________________
211 AliEMCALDigitizer::~AliEMCALDigitizer()
212{
14ce0a6e 213 //dtor
5dee926e 214 if (AliRunLoader::GetRunLoader()) {
215 AliLoader *emcalLoader=0;
216 if ((emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")))
217 emcalLoader->CleanDigitizer();
218 }
219 else
220 AliDebug(1," no runloader present");
f5f384f4 221 delete [] fInputFileNames ;
222 delete [] fEventNames ;
64942713 223
61e0abb5 224}
225
226//____________________________________________________________________________
09884213 227void AliEMCALDigitizer::Digitize(Int_t event)
05a92d59 228{
61e0abb5 229
230 // Makes the digitization of the collected summable digits
231 // for this it first creates the array of all EMCAL modules
ab6a174f 232 // filled with noise and after that adds contributions from
233 // SDigits. This design helps to avoid scanning over the
234 // list of digits to add contribution of any new SDigit.
235 //
236 // JLK 26-Jun-2008
237 // Note that SDigit energy info is stored as an amplitude, so we
238 // must call the Calibrate() method of the SDigitizer to convert it
239 // back to an energy in GeV before adding it to the Digit
240 //
1963b290 241 static int nEMC=0; //max number of digits possible
61e0abb5 242
5dee926e 243 AliRunLoader *rl = AliRunLoader::GetRunLoader();
244 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
14ce0a6e 245 Int_t readEvent = event ;
1963b290 246 // fManager is data member from AliDigitizer
45fa49ca 247 if (fManager)
14ce0a6e 248 readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
cf24bdc4 249 AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
14ce0a6e 250 readEvent, GetTitle(), fEventFolderName.Data())) ;
251 rl->GetEvent(readEvent);
5dee926e 252
253 TClonesArray * digits = emcalLoader->Digits() ;
ab6a174f 254 digits->Delete() ; //correct way to clear array when memory is
255 //allocated by objects stored in array
61e0abb5 256
e4dcc484 257 // Load Geometry
30bc5594 258 AliEMCALGeometry *geom = 0;
259 if (rl->GetAliRun()) {
260 AliEMCAL * emcal = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
261 geom = emcal->GetGeometry();
262 }
263 else
264 AliFatal("Could not get AliRun from runLoader");
e4dcc484 265
fe93d365 266 nEMC = geom->GetNCells();
267 AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
268
88cb7938 269 Int_t absID ;
61e0abb5 270
88cb7938 271 digits->Expand(nEMC) ;
05a92d59 272
88cb7938 273 // get first the sdigitizer from the tasks list (must have same name as the digitizer)
5dee926e 274 if ( !emcalLoader->SDigitizer() )
275 emcalLoader->LoadSDigitizer();
276 AliEMCALSDigitizer * sDigitizer = dynamic_cast<AliEMCALSDigitizer*>(emcalLoader->SDigitizer());
88cb7938 277
278 if ( !sDigitizer )
279 Fatal("Digitize", "SDigitizer with name %s %s not found", fEventFolderName.Data(), GetTitle() ) ;
280
281 //take all the inputs to add together and load the SDigits
282 TObjArray * sdigArray = new TObjArray(fInput) ;
5dee926e 283 sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
88cb7938 284 Int_t i ;
5dee926e 285
88cb7938 286 for(i = 1 ; i < fInput ; i++){
287 TString tempo(fEventNames[i]) ;
288 tempo += i ;
3880b0fb 289
290 AliRunLoader * rl2 = AliRunLoader::GetRunLoader(tempo) ;
291
292 if (rl2==0)
293 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
294
45fa49ca 295 if (fManager)
14ce0a6e 296 readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
297 Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
3880b0fb 298 rl2->LoadSDigits();
14ce0a6e 299 rl2->GetEvent(readEvent);
3880b0fb 300 AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
301 sdigArray->AddAt(emcalLoader2->SDigits(), i) ;
0c219753 302 }
5dee926e 303
d968cee0 304 //Find the first tower with signal
a9485c37 305 Int_t nextSig = nEMC + 1 ;
88cb7938 306 TClonesArray * sdigits ;
307 for(i = 0 ; i < fInput ; i++){
308 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
814ad4bf 309 if ( !sdigits->GetEntriesFast() )
310 continue ;
88cb7938 311 Int_t curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(0))->GetId() ;
814ad4bf 312 if(curNext < nextSig)
313 nextSig = curNext ;
cf24bdc4 314 AliDebug(1,Form("input %i : #sdigits %i \n",
315 i, sdigits->GetEntriesFast()));
814ad4bf 316 }
cf24bdc4 317 AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
61e0abb5 318
88cb7938 319 TArrayI index(fInput) ;
814ad4bf 320 index.Reset() ; //Set all indexes to zero
61e0abb5 321
05a92d59 322 AliEMCALDigit * digit ;
323 AliEMCALDigit * curSDigit ;
ffa6d63b 324
9517d886 325 // TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ;
ffa6d63b 326
5dee926e 327 //Put Noise contribution
d25f2c54 328 for(absID = 0; absID < nEMC; absID++){ // Nov 30, 2006 by PAI; was from 1 to nEMC
ab6a174f 329 Float_t energy = 0 ;
a9485c37 330 // amplitude set to zero, noise will be added later
d25f2c54 331 new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0, TimeOfNoise() ); // absID-1->absID
814ad4bf 332 //look if we have to add signal?
d25f2c54 333 digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
88cb7938 334
814ad4bf 335 if(absID==nextSig){
54cbf5a6 336 //Add SDigits from all inputs
9517d886 337 // ticks->Clear() ;
338 //Int_t contrib = 0 ;
339
340 //Follow PHOS and comment out this timing model til a better one
341 //can be developed - JLK 28-Apr-2008
342
343 //Float_t a = digit->GetAmp() ;
344 //Float_t b = TMath::Abs( a /fTimeSignalLength) ;
88cb7938 345 //Mark the beginning of the signal
9517d886 346 //new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime(),0, b);
fdebddeb 347 //Mark the end of the signal
9517d886 348 //new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime()+fTimeSignalLength, -a, -b);
349
350 // Calculate time as time of the largest digit
351 Float_t time = digit->GetTime() ;
ab6a174f 352 Float_t aTime= digit->GetAmp() ;
814ad4bf 353
88cb7938 354 // loop over input
355 for(i = 0; i< fInput ; i++){ //loop over (possible) merge sources
356 if(dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
357 curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
814ad4bf 358 else
359 curSDigit = 0 ;
360 //May be several digits will contribute from the same input
fdebddeb 361 while(curSDigit && (curSDigit->GetId() == absID)){
814ad4bf 362 //Shift primary to separate primaries belonging different inputs
363 Int_t primaryoffset ;
364 if(fManager)
365 primaryoffset = fManager->GetMask(i) ;
366 else
367 primaryoffset = i ;
368 curSDigit->ShiftPrimary(primaryoffset) ;
9517d886 369
370 //Remove old timing model - JLK 28-April-2008
371 //a = curSDigit->GetAmp() ;
372 //b = a /fTimeSignalLength ;
373 //new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime(),0, b);
374 //new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
ab6a174f 375 if(curSDigit->GetAmp()>aTime) {
376 aTime = curSDigit->GetAmp();
9517d886 377 time = curSDigit->GetTime();
378 }
814ad4bf 379
ab6a174f 380 *digit = *digit + *curSDigit ; //adds amplitudes of each digit
814ad4bf 381
382 index[i]++ ;
88cb7938 383 if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
1963b290 384 curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
814ad4bf 385 else
386 curSDigit = 0 ;
387 }
388 }
ab6a174f 389 //Here we convert the summed amplitude to an energy in GeV
390 energy = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
a9485c37 391 // add fluctuations for photo-electron creation
ab6a174f 392 energy *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
88cb7938 393
05a92d59 394 //calculate and set time
9517d886 395 //New timing model needed - JLK 28-April-2008
396 //Float_t time = FrontEdgeTime(ticks) ;
814ad4bf 397 digit->SetTime(time) ;
398
399 //Find next signal module
a9485c37 400 nextSig = nEMC + 1 ;
88cb7938 401 for(i = 0 ; i < fInput ; i++){
402 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
814ad4bf 403 Int_t curNext = nextSig ;
404 if(sdigits->GetEntriesFast() > index[i] ){
88cb7938 405 curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]))->GetId() ;
814ad4bf 406 }
407 if(curNext < nextSig) nextSig = curNext ;
408 }
409 }
a9485c37 410 // add the noise now
ab6a174f 411 energy += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
412 // JLK 26-June-2008
413 //Now digitize the energy according to the sDigitizer method,
414 //which merely converts the energy to an integer. Later we will
415 //check that the stored value matches our allowed dynamic ranges
416 digit->SetAmp(sDigitizer->Digitize(energy)) ;
417 AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
418 absID, energy, nextSig));
419 } // for(absID = 0; absID < nEMC; absID++)
814ad4bf 420
9517d886 421 //ticks->Delete() ;
422 //delete ticks ;
61e0abb5 423
ab6a174f 424 //JLK is it better to call Clear() here?
05a92d59 425 delete sdigArray ; //We should not delete its contents
61e0abb5 426
814ad4bf 427 //remove digits below thresholds
88cb7938 428 for(i = 0 ; i < nEMC ; i++){
429 digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
ab6a174f 430 Float_t threshold = fDigitThreshold ; //this is in GeV
431 //need to calibrate digit amplitude to energy in GeV for comparison
88cb7938 432 if(sDigitizer->Calibrate( digit->GetAmp() ) < threshold)
433 digits->RemoveAt(i) ;
434 else
814ad4bf 435 digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
a9485c37 436 }
61e0abb5 437
814ad4bf 438 digits->Compress() ;
61e0abb5 439
05a92d59 440 Int_t ndigits = digits->GetEntriesFast() ;
7ea6391b 441
ab6a174f 442 //JLK 26-June-2008
443 //After we have done the summing and digitizing to create the
444 //digits, now we want to calibrate the resulting amplitude to match
445 //the dynamic range of our real data.
446 Float_t energy=0;
05a92d59 447 for (i = 0 ; i < ndigits ; i++) {
448 digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
814ad4bf 449 digit->SetIndexInList(i) ;
1963b290 450 energy = sDigitizer->Calibrate(digit->GetAmp()) ;
ab6a174f 451 digit->SetAmp(DigitizeEnergy(energy, digit->GetId()) ) ;
61e0abb5 452 }
ab6a174f 453
814ad4bf 454}
61e0abb5 455
d62a891a 456// //_____________________________________________________________________
e4dcc484 457Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
814ad4bf 458{
ab6a174f 459 // JLK 26-June-2008
e4dcc484 460 // Returns digitized value of the energy in a cell absId
ab6a174f 461 // using the calibration constants stored in the OCDB
462 // or default values if no CalibData object is found.
463 // This effectively converts everything to match the dynamic range
464 // of the real data we will collect
465 //
53f1c378 466 // Load Geometry
467 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
e4dcc484 468
469 if (geom==0)
d25f2c54 470 AliFatal("Did not get geometry from EMCALLoader");
e4dcc484 471
472 Int_t iSupMod = -1;
2bb3725c 473 Int_t nModule = -1;
e4dcc484 474 Int_t nIphi = -1;
475 Int_t nIeta = -1;
476 Int_t iphi = -1;
477 Int_t ieta = -1;
fdebddeb 478 Int_t channel = -999;
e4dcc484 479
2bb3725c 480 Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
e4dcc484 481 if(!bCell)
d25f2c54 482 Error("DigitizeEnergy","Wrong cell id number : AbsId %i ", AbsId) ;
2bb3725c 483 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
d62a891a 484
53f1c378 485 if(fCalibData) {
486 fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
487 fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
e4dcc484 488 }
5eb74a51 489
490 channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC )) ;
491
fdebddeb 492 if(channel > fNADCEC )
7d74eaca 493 channel = fNADCEC ;
814ad4bf 494 return channel ;
d62a891a 495
61e0abb5 496}
497
498//____________________________________________________________________________
05a92d59 499void AliEMCALDigitizer::Exec(Option_t *option)
500{
45fa49ca 501 // Steering method to process digitization for events
502 // in the range from fFirstEvent to fLastEvent.
503 // This range is optionally set by SetEventRange().
504 // if fLastEvent=-1, then process events until the end.
505 // by default fLastEvent = fFirstEvent (process only one event)
05a92d59 506
88cb7938 507 if (!fInit) { // to prevent overwrite existing file
508 Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
509 return ;
510 }
511
814ad4bf 512 if (strstr(option,"print")) {
b42d4197 513
88cb7938 514 Print();
814ad4bf 515 return ;
516 }
517
61e0abb5 518 if(strstr(option,"tim"))
519 gBenchmark->Start("EMCALDigitizer");
5dee926e 520
521 AliRunLoader *rl = AliRunLoader::GetRunLoader();
522 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
523
5b4d2c50 524 // Post Digitizer to the white board
5dee926e 525 emcalLoader->PostDigitizer(this) ;
5b4d2c50 526
30bc5594 527 if (fLastEvent == -1) {
5dee926e 528 fLastEvent = rl->GetNumberOfEvents() - 1 ;
30bc5594 529 }
396a348e 530 else if (fManager)
1963b290 531 fLastEvent = fFirstEvent ; // what is this ??
396a348e 532
45fa49ca 533 Int_t nEvents = fLastEvent - fFirstEvent + 1;
cf24bdc4 534 Int_t ievent;
535
5dee926e 536 rl->LoadSDigits("EMCAL");
cf24bdc4 537 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
64942713 538
5dee926e 539 rl->GetEvent(ievent);
05a92d59 540
814ad4bf 541 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
88cb7938 542
9e5d2067 543 WriteDigits() ;
88cb7938 544
61e0abb5 545 if(strstr(option,"deb"))
546 PrintDigits(option);
1963b290 547 if(strstr(option,"table")) gObjectTable->Print();
88cb7938 548
814ad4bf 549 //increment the total number of Digits per run
5dee926e 550 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
61e0abb5 551 }
64942713 552
5dee926e 553 emcalLoader->CleanDigitizer() ;
5b4d2c50 554
61e0abb5 555 if(strstr(option,"tim")){
556 gBenchmark->Stop("EMCALDigitizer");
b42d4197 557 AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event",
558 gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
88cb7938 559 }
61e0abb5 560}
561
05a92d59 562//____________________________________________________________________________
563Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks)
88cb7938 564{
565 // Returns the shortest time among all time ticks
566
05a92d59 567 ticks->Sort() ; //Sort in accordance with times of ticks
568 TIter it(ticks) ;
569 AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
570 Float_t time = ctick->CrossingTime(fTimeThreshold) ;
571
572 AliEMCALTick * t ;
573 while((t=(AliEMCALTick*) it.Next())){
574 if(t->GetTime() < time) //This tick starts before crossing
575 *ctick+=*t ;
576 else
577 return time ;
578
579 time = ctick->CrossingTime(fTimeThreshold) ;
580 }
581 return time ;
582}
583
584//____________________________________________________________________________
585Bool_t AliEMCALDigitizer::Init()
586{
587 // Makes all memory allocations
88cb7938 588 fInit = kTRUE ;
5dee926e 589 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
590
591 if ( emcalLoader == 0 ) {
592 Fatal("Init", "Could not obtain the AliEMCALLoader");
05a92d59 593 return kFALSE;
594 }
05a92d59 595
45fa49ca 596 fFirstEvent = 0 ;
597 fLastEvent = fFirstEvent ;
598
88cb7938 599 if(fManager)
600 fInput = fManager->GetNinputs() ;
601 else
602 fInput = 1 ;
603
604 fInputFileNames = new TString[fInput] ;
605 fEventNames = new TString[fInput] ;
606 fInputFileNames[0] = GetTitle() ;
607 fEventNames[0] = fEventFolderName.Data() ;
608 Int_t index ;
609 for (index = 1 ; index < fInput ; index++) {
610 fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
611 TString tempo = fManager->GetInputFolderName(index) ;
612 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager
05a92d59 613 }
614
88cb7938 615 //to prevent cleaning of this object while GetEvent is called
5dee926e 616 emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
e52475ed 617
53f1c378 618 //Calibration instance
619 fCalibData = emcalLoader->CalibData();
88cb7938 620 return fInit ;
05a92d59 621}
622
623//____________________________________________________________________________
624void AliEMCALDigitizer::InitParameters()
14ce0a6e 625{
29b5d492 626 // Parameter initialization for digitizer
627 // Tune parameters - 24-nov-04; Apr 29, 2007
fe93d365 628 // New parameters JLK 14-Apr-2008
1963b290 629
fe93d365 630 fMeanPhotonElectron = 4400; // electrons per GeV
631 fPinNoise = 0.037; // pin noise in GEV from analysis test beam data
88cb7938 632 if (fPinNoise == 0. )
633 Warning("InitParameters", "No noise added\n") ;
1963b290 634 fDigitThreshold = fPinNoise * 3; // 3 * sigma
635 fTimeResolution = 0.3e-9 ; // 300 psc
05a92d59 636 fTimeSignalLength = 1.0e-9 ;
05a92d59 637
97075cd8 638 // These defaults are normally not used.
639 // Values are read from calibration database instead
640 fADCchannelEC = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
641 fADCpedestalEC = 0.0 ; // GeV
1963b290 642 fNADCEC = (Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
a9485c37 643
1963b290 644 fTimeThreshold = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
ab6a174f 645
05a92d59 646}
61e0abb5 647
61e0abb5 648//__________________________________________________________________
09884213 649void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
05a92d59 650{
651 // Allows to produce digits by superimposing background and signal event.
61e0abb5 652 // It is assumed, that headers file with SIGNAL events is opened in
05a92d59 653 // the constructor.
654 // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed
655 // Thus we avoid writing (changing) huge and expensive
61e0abb5 656 // backgound files: all output will be writen into SIGNAL, i.e.
657 // opened in constructor file.
658 //
05a92d59 659 // One can open as many files to mix with as one needs.
660 // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
661 // can be mixed.
61e0abb5 662
05a92d59 663 if( strcmp(GetName(), "") == 0 )
61e0abb5 664 Init() ;
814ad4bf 665
05a92d59 666 if(fManager){
9859bfc0 667 Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
61e0abb5 668 return ;
05a92d59 669 }
88cb7938 670 // looking for file which contains AliRun
671 if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
672 Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
673 return ;
61e0abb5 674 }
88cb7938 675 // looking for the file which contains SDigits
5dee926e 676 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
677 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
e191bb57 678 if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
88cb7938 679 fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
680 if ( (gSystem->AccessPathName(fileName)) ) {
681 Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
682 return ;
814ad4bf 683 }
88cb7938 684 // need to increase the arrays
30bc5594 685 // MvL: This code only works when fInput == 1, I think.
88cb7938 686 TString tempo = fInputFileNames[fInput-1] ;
687 delete [] fInputFileNames ;
688 fInputFileNames = new TString[fInput+1] ;
689 fInputFileNames[fInput-1] = tempo ;
690
691 tempo = fEventNames[fInput-1] ;
692 delete [] fEventNames ;
693 fEventNames = new TString[fInput+1] ;
694 fEventNames[fInput-1] = tempo ;
695
696 fInputFileNames[fInput] = alirunFileName ;
697 fEventNames[fInput] = eventFolderName ;
698 fInput++ ;
699}
1963b290 700
fe93d365 701//__________________________________________________________________
1963b290 702void AliEMCALDigitizer::Print1(Option_t * option)
703{ // 19-nov-04 - just for convinience
704 Print();
705 PrintDigits(option);
706}
707
61e0abb5 708//__________________________________________________________________
eb0b1051 709void AliEMCALDigitizer::Print(Option_t*)const
88cb7938 710{
711 // Print Digitizer's parameters
fdebddeb 712 printf("Print: \n------------------- %s -------------", GetName() ) ;
88cb7938 713 if( strcmp(fEventFolderName.Data(), "") != 0 ){
714 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
715
716 Int_t nStreams ;
717 if (fManager)
718 nStreams = GetNInputStreams() ;
719 else
720 nStreams = fInput ;
721
722 Int_t index = 0 ;
5dee926e 723
724 AliRunLoader *rl=0;
725
88cb7938 726 for (index = 0 ; index < nStreams ; index++) {
727 TString tempo(fEventNames[index]) ;
728 tempo += index ;
5dee926e 729 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
730 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
731 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
732 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
e191bb57 733 if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
88cb7938 734 fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
735 printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
736 }
5dee926e 737
738 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
739
740 printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
88cb7938 741
742 printf("\nWith following parameters:\n") ;
743
744 printf(" Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise) ;
745 printf(" Threshold in EMC (fDigitThreshold) = %f\n", fDigitThreshold) ;
746 printf("---------------------------------------------------\n") ;
ea5f3389 747 }
748 else
fdebddeb 749 printf("Print: AliEMCALDigitizer not initialized") ;
61e0abb5 750}
814ad4bf 751
61e0abb5 752//__________________________________________________________________
14ce0a6e 753void AliEMCALDigitizer::PrintDigits(Option_t * option)
754{
755 //utility method for printing digit information
756
5dee926e 757 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
758 TClonesArray * digits = emcalLoader->Digits() ;
759 TClonesArray * sdigits = emcalLoader->SDigits() ;
88cb7938 760
1963b290 761 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
5dee926e 762 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
88cb7938 763
764 if(strstr(option,"all")){
61e0abb5 765 //loop over digits
766 AliEMCALDigit * digit;
88cb7938 767 printf("\nEMC digits (with primaries):\n") ;
768 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
61e0abb5 769 Int_t index ;
88cb7938 770 for (index = 0 ; index < digits->GetEntries() ; index++) {
771 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
772 printf("\n%6d %8d %6.5e %4d %2d : ",
11f9c5ff 773 digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
61e0abb5 774 Int_t iprimary;
9859bfc0 775 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
88cb7938 776 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
9859bfc0 777 }
778 }
61e0abb5 779 }
1963b290 780 printf("\n");
61e0abb5 781}
05a92d59 782
783//__________________________________________________________________
784Float_t AliEMCALDigitizer::TimeOfNoise(void)
04f0bda3 785{
786 // Calculates the time signal generated by noise
26a2ef9d 787 //PH Info("TimeOfNoise", "Change me") ;
04f0bda3 788 return gRandom->Rndm() * 1.28E-5;
05a92d59 789}
790
88cb7938 791//__________________________________________________________________
792void AliEMCALDigitizer::Unload()
793{
fdebddeb 794 // Unloads the SDigits and Digits
5dee926e 795 AliRunLoader *rl=0;
796
88cb7938 797 Int_t i ;
798 for(i = 1 ; i < fInput ; i++){
799 TString tempo(fEventNames[i]) ;
5dee926e 800 tempo += i;
801 if ((rl = AliRunLoader::GetRunLoader(tempo)))
802 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
88cb7938 803 }
5dee926e 804 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
805 emcalLoader->UnloadDigits() ;
88cb7938 806}
807
814ad4bf 808//_________________________________________________________________________________________
9e5d2067 809void AliEMCALDigitizer::WriteDigits()
814ad4bf 810{
61e0abb5 811
814ad4bf 812 // Makes TreeD in the output file.
813 // Check if branch already exists:
814 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
815 // already existing branches.
816 // else creates branch with Digits, named "EMCAL", title "...",
817 // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
818 // and names of files, from which digits are made.
61e0abb5 819
5dee926e 820 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
cf24bdc4 821
5dee926e 822 const TClonesArray * digits = emcalLoader->Digits() ;
823 TTree * treeD = emcalLoader->TreeD();
824 if ( !treeD ) {
825 emcalLoader->MakeDigitsContainer();
826 treeD = emcalLoader->TreeD();
827 }
88cb7938 828
829 // -- create Digits branch
830 Int_t bufferSize = 32000 ;
5dee926e 831 TBranch * digitsBranch = 0;
30bc5594 832 if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
5dee926e 833 digitsBranch->SetAddress(&digits);
30bc5594 834 AliWarning("Digits Branch already exists. Not all digits will be visible");
835 }
5dee926e 836 else
837 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
838 //digitsBranch->SetTitle(fEventFolderName);
839 treeD->Fill() ;
88cb7938 840
5dee926e 841 emcalLoader->WriteDigits("OVERWRITE");
842 emcalLoader->WriteDigitizer("OVERWRITE");
88cb7938 843
844 Unload() ;
845
814ad4bf 846}
61e0abb5 847
fe93d365 848//__________________________________________________________________
1963b290 849void AliEMCALDigitizer::Browse(TBrowser* b)
850{
1963b290 851 TTask::Browse(b);
852}