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