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