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