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