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