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