new file to test SDigits
[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//
6cc75819 23// In addition it performs mixing/embedding of summable digits from different events.
61e0abb5 24//
61d80ce0 25// For each event 3 branches are created in TreeD:
61e0abb5 26// "EMCAL" - list of digits
61d80ce0 27// "EMCALTRG" - list of trigger digits
61e0abb5 28// "AliEMCALDigitizer" - AliEMCALDigitizer with all parameters used in digitization
29//
61e0abb5 30//
ffa6d63b 31////////////////////////////////////////////////////////////////////////////////////
61e0abb5 32//
ffa6d63b 33//*-- Author: Sahal Yacoob (LBL)
814ad4bf 34// based on : AliEMCALDigitizer
05a92d59 35// Modif:
36// August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
6cc75819 37// of new IO (a la PHOS)
1963b290 38// November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry
61d80ce0 39// July 2011 GCB: Digitizer modified to accomodate embedding.
40// Time calibration added. Decalibration possibility of energy and time added
1963b290 41//_________________________________________________________________________________
61e0abb5 42
43// --- ROOT system ---
e939a978 44#include <TROOT.h>
45#include <TTree.h>
46#include <TSystem.h>
47#include <TBenchmark.h>
e939a978 48#include <TBrowser.h>
49#include <TObjectTable.h>
50#include <TRandom.h>
916f1e76 51#include <TF1.h>
21e7df44 52#include <cassert>
61e0abb5 53
54// --- AliRoot header files ---
b42d4197 55#include "AliLog.h"
e4dcc484 56#include "AliRun.h"
f21fc003 57#include "AliDigitizationInput.h"
5dee926e 58#include "AliRunLoader.h"
53f1c378 59#include "AliCDBManager.h"
60#include "AliCDBEntry.h"
61e0abb5 61#include "AliEMCALDigit.h"
05a92d59 62#include "AliEMCAL.h"
5dee926e 63#include "AliEMCALLoader.h"
61e0abb5 64#include "AliEMCALDigitizer.h"
65#include "AliEMCALSDigitizer.h"
66#include "AliEMCALGeometry.h"
a435f763 67#include "AliEMCALCalibData.h"
33a52e55 68#include "AliEMCALSimParam.h"
916f1e76 69#include "AliEMCALRawDigit.h"
6995e8b7 70#include "AliCaloCalibPedestal.h"
916f1e76 71
61d80ce0 72 namespace
73 {
74 Double_t HeavisideTheta(Double_t x)
75 {
76 Double_t signal = 0.;
77
78 if (x > 0.) signal = 1.;
79
80 return signal;
81 }
82
83 Double_t AnalogFastORFunction(Double_t *x, Double_t *par)
84 {
85 Double_t v0 = par[0];
86 Double_t t0 = par[1];
87 Double_t tr = par[2];
88
89 Double_t R1 = 1000.;
90 Double_t C1 = 33e-12;
91 Double_t R2 = 1800;
92 Double_t C2 = 22e-12;
93
94 Double_t t = x[0];
95
96 return (((0.8*(-((TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-t + t0)/(C1*R1))*
97 TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2)) +
98 C1*C2*R1*R2*(1 - (C2*TMath::Power(TMath::E(),(-t + t0)/(C2*R2))*R2)/(-(C1*R1) + C2*R2)))*v0*
99 HeavisideTheta(t - t0))/tr
100 - (0.8*(C1*C2*R1*R2 -
101 (TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C1*R1))*
102 TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2) +
103 (C1*TMath::Power(C2,2)*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C2*R2))*
104 R1*TMath::Power(R2,2))/(C1*R1 - C2*R2))*v0*
105 HeavisideTheta(t - t0 - 1.25*tr))/tr)/(C2*R1));
106 }
107 }
05a92d59 108
61e0abb5 109ClassImp(AliEMCALDigitizer)
61d80ce0 110
111
61e0abb5 112//____________________________________________________________________________
18a21c7c 113AliEMCALDigitizer::AliEMCALDigitizer()
114 : AliDigitizer("",""),
115 fDefaultInit(kTRUE),
116 fDigitsInRun(0),
117 fInit(0),
118 fInput(0),
119 fInputFileNames(0x0),
120 fEventNames(0x0),
121 fDigitThreshold(0),
122 fMeanPhotonElectron(0),
5fb45456 123 fGainFluctuations(0),
18a21c7c 124 fPinNoise(0),
6cc75819 125 fTimeNoise(0),
f0a6dc6f 126 fTimeDelay(0),
a2f8e711 127 fTimeResolutionPar0(0),
128 fTimeResolutionPar1(0),
18a21c7c 129 fADCchannelEC(0),
130 fADCpedestalEC(0),
6cc75819 131 fADCchannelECDecal(0),
132 fTimeChannel(0),
133 fTimeChannelDecal(0),
18a21c7c 134 fNADCEC(0),
135 fEventFolderName(""),
136 fFirstEvent(0),
137 fLastEvent(0),
f21fc003 138 fCalibData(0x0),
139 fSDigitizer(0x0)
61e0abb5 140{
141 // ctor
88cb7938 142 InitParameters() ;
f21fc003 143 fDigInput = 0 ; // We work in the standalone mode
839828a6 144}
61e0abb5 145
839828a6 146//____________________________________________________________________________
18a21c7c 147AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
f21fc003 148 : AliDigitizer("EMCALDigitizer", alirunFileName),
18a21c7c 149 fDefaultInit(kFALSE),
150 fDigitsInRun(0),
151 fInit(0),
152 fInput(0),
153 fInputFileNames(0),
154 fEventNames(0),
155 fDigitThreshold(0),
156 fMeanPhotonElectron(0),
5fb45456 157 fGainFluctuations(0),
18a21c7c 158 fPinNoise(0),
6cc75819 159 fTimeNoise(0),
96d5ea0d 160 fTimeDelay(0),
a2f8e711 161 fTimeResolutionPar0(0),
162 fTimeResolutionPar1(0),
18a21c7c 163 fADCchannelEC(0),
164 fADCpedestalEC(0),
6cc75819 165 fADCchannelECDecal(0),
166 fTimeChannel(0),
167 fTimeChannelDecal(0),
18a21c7c 168 fNADCEC(0),
169 fEventFolderName(eventFolderName),
170 fFirstEvent(0),
171 fLastEvent(0),
f21fc003 172 fCalibData(0x0),
173 fSDigitizer(0x0)
839828a6 174{
05a92d59 175 // ctor
a9485c37 176 InitParameters() ;
88cb7938 177 Init() ;
f21fc003 178 fDigInput = 0 ; // We work in the standalone mode
61e0abb5 179}
839828a6 180
61e0abb5 181//____________________________________________________________________________
18a21c7c 182AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
183 : AliDigitizer(d.GetName(),d.GetTitle()),
184 fDefaultInit(d.fDefaultInit),
185 fDigitsInRun(d.fDigitsInRun),
186 fInit(d.fInit),
187 fInput(d.fInput),
188 fInputFileNames(d.fInputFileNames),
189 fEventNames(d.fEventNames),
190 fDigitThreshold(d.fDigitThreshold),
191 fMeanPhotonElectron(d.fMeanPhotonElectron),
5fb45456 192 fGainFluctuations(d.fGainFluctuations),
18a21c7c 193 fPinNoise(d.fPinNoise),
6cc75819 194 fTimeNoise(d.fTimeNoise),
f0a6dc6f 195 fTimeDelay(d.fTimeDelay),
a2f8e711 196 fTimeResolutionPar0(d.fTimeResolutionPar0),
197 fTimeResolutionPar1(d.fTimeResolutionPar1),
18a21c7c 198 fADCchannelEC(d.fADCchannelEC),
199 fADCpedestalEC(d.fADCpedestalEC),
6cc75819 200 fADCchannelECDecal(d.fADCchannelECDecal),
201 fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
18a21c7c 202 fNADCEC(d.fNADCEC),
203 fEventFolderName(d.fEventFolderName),
204 fFirstEvent(d.fFirstEvent),
205 fLastEvent(d.fLastEvent),
f21fc003 206 fCalibData(d.fCalibData),
207 fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
88cb7938 208{
209 // copyy ctor
96d5ea0d 210}
88cb7938 211
212//____________________________________________________________________________
f21fc003 213AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
214 : AliDigitizer(rd,"EMCALDigitizer"),
18a21c7c 215 fDefaultInit(kFALSE),
216 fDigitsInRun(0),
217 fInit(0),
218 fInput(0),
219 fInputFileNames(0),
220 fEventNames(0),
da509d54 221 fDigitThreshold(0),
18a21c7c 222 fMeanPhotonElectron(0),
5fb45456 223 fGainFluctuations(0),
da509d54 224 fPinNoise(0.),
6cc75819 225 fTimeNoise(0.),
f0a6dc6f 226 fTimeDelay(0.),
a2f8e711 227 fTimeResolutionPar0(0.),
228 fTimeResolutionPar1(0.),
18a21c7c 229 fADCchannelEC(0),
230 fADCpedestalEC(0),
6cc75819 231 fADCchannelECDecal(0),
232 fTimeChannel(0),
233 fTimeChannelDecal(0),
18a21c7c 234 fNADCEC(0),
235 fEventFolderName(0),
236 fFirstEvent(0),
237 fLastEvent(0),
f21fc003 238 fCalibData(0x0),
239 fSDigitizer(0x0)
61e0abb5 240{
45fa49ca 241 // ctor Init() is called by RunDigitizer
f21fc003 242 fDigInput = rd ;
243 fEventFolderName = fDigInput->GetInputFolderName(0) ;
244 SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
05a92d59 245 InitParameters() ;
839828a6 246}
814ad4bf 247
839828a6 248//____________________________________________________________________________
249 AliEMCALDigitizer::~AliEMCALDigitizer()
250{
14ce0a6e 251 //dtor
f5f384f4 252 delete [] fInputFileNames ;
253 delete [] fEventNames ;
f21fc003 254 if (fSDigitizer) delete fSDigitizer;
61e0abb5 255}
256
257//____________________________________________________________________________
a81281f7 258void AliEMCALDigitizer::Digitize(Int_t event)
259{
61e0abb5 260 // Makes the digitization of the collected summable digits
261 // for this it first creates the array of all EMCAL modules
a81281f7 262 // filled with noise and after that adds contributions from
263 // SDigits. This design helps to avoid scanning over the
ab6a174f 264 // list of digits to add contribution of any new SDigit.
265 //
266 // JLK 26-Jun-2008
267 // Note that SDigit energy info is stored as an amplitude, so we
268 // must call the Calibrate() method of the SDigitizer to convert it
269 // back to an energy in GeV before adding it to the Digit
270 //
1963b290 271 static int nEMC=0; //max number of digits possible
33c3c91a 272 AliRunLoader *rl = AliRunLoader::Instance();
5dee926e 273 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
7e1d9a9b 274
a81281f7 275 if(emcalLoader)
276 {
277 Int_t readEvent = event ;
f21fc003 278 // fDigInput is data member from AliDigitizer
a81281f7 279 if (fDigInput)
280 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
281 AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
282 readEvent, GetTitle(), fEventFolderName.Data())) ;
6cc75819 283 rl->GetEvent(readEvent);
7e1d9a9b 284
a81281f7 285 TClonesArray * digits = emcalLoader->Digits() ;
6cc75819 286 digits->Delete() ; //correct way to clear array when memory is
287 //allocated by objects stored in array
88cb7938 288
6cc75819 289 // Load Geometry
290 AliEMCALGeometry *geom = 0;
a81281f7 291 if (!rl->GetAliRun())
292 {
6cc75819 293 AliFatal("Could not get AliRun from runLoader");
96d5ea0d 294 }
6cc75819 295 else{
296 AliEMCAL * emcal = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
297 geom = emcal->GetGeometry();
298 nEMC = geom->GetNCells();
299 AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
300 }
301
302 Int_t absID = -1 ;
303
304 digits->Expand(nEMC) ;
96d5ea0d 305
f21fc003 306 // RS create a digitizer on the fly
307 if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
308 fSDigitizer->SetEventRange(0, -1) ;
309 //
310 //take all the inputs to add together and load the SDigits
311 TObjArray * sdigArray = new TObjArray(fInput) ;
312 sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
313 Int_t i ;
314 Int_t embed = kFALSE;
a81281f7 315 for(i = 1 ; i < fInput ; i++)
316 {
317 TString tempo(fEventNames[i]) ;
f21fc003 318 tempo += i ;
a81281f7 319
320 AliRunLoader * rl2 = AliRunLoader::GetRunLoader(tempo) ;
321
322 if (rl2==0)
323 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
324
325 if (fDigInput)
326 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
327 Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
f21fc003 328 rl2->LoadSDigits();
329 // rl2->LoadDigits();
330 rl2->GetEvent(readEvent);
331 if(rl2->GetDetectorLoader("EMCAL"))
a81281f7 332 {
333 AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
334
335 if(emcalLoader2)
6cc75819 336 {
a81281f7 337 //Merge 2 simulated sdigits
338 if(emcalLoader2->SDigits()){
339 TClonesArray* sdigits2 = emcalLoader2->SDigits();
340 sdigArray->AddAt(sdigits2, i) ;
341 // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
342 // do not smear energy of embedded, do not add noise to any sdigits
343 if(sdigits2->GetEntriesFast()>0)
344 {
345 //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
346 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
347 if( digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded ) embed = kTRUE;
6cc75819 348 }
a81281f7 349 }
6cc75819 350
a81281f7 351 }//loader2
352 else AliFatal("EMCAL Loader of second event not available!");
353
354 }
f21fc003 355 else Fatal("Digitize", "Loader of second input not found");
356 }// input loop
a81281f7 357
f21fc003 358 //Find the first tower with signal
a81281f7 359 Int_t nextSig = nEMC + 1 ;
360 TClonesArray * sdigits ;
f21fc003 361 for(i = 0 ; i < fInput ; i++){
362 if(i > 0 && embed) continue;
363 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
a81281f7 364 if(sdigits)
365 {
366 if (sdigits->GetEntriesFast() )
367 {
368 AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
369 if(sd)
370 {
371 Int_t curNext = sd->GetId() ;
372 if(curNext < nextSig)
373 nextSig = curNext ;
374 AliDebug(1,Form("input %i : #sdigits %i \n",
375 i, sdigits->GetEntriesFast()));
376
377 }//First entry is not NULL
378 else
379 {
380 AliDebug(1, "NULL sdigit pointer");
381 continue;
382 }
383 }//There is at least one entry
384 else
385 {
386 AliDebug(1, "NULL sdigits array");
387 continue;
388 }
f21fc003 389 }// SDigits array exists
a81281f7 390 else
391 {
392 AliDebug(1,"Null sdigit pointer");
393 continue;
f21fc003 394 }
395 }// input loop
396 AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
a81281f7 397
f21fc003 398 TArrayI index(fInput) ;
399 index.Reset() ; //Set all indexes to zero
a81281f7 400
f21fc003 401 AliEMCALDigit * digit ;
402 AliEMCALDigit * curSDigit ;
a81281f7 403
404 //Put Noise contribution, smear time and energy
f21fc003 405 Float_t timeResolution = 0;
a81281f7 406 for(absID = 0; absID < nEMC; absID++)
407 { // Nov 30, 2006 by PAI; was from 1 to nEMC
408
6995e8b7 409 if (IsDead(absID)) continue; // Don't instantiate dead digits
a81281f7 410
411 Float_t energy = 0 ;
412
f21fc003 413 // amplitude set to zero, noise will be added later
414 Float_t noiseTime = 0.;
415 if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
416 new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
417 //look if we have to add signal?
418 digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
a81281f7 419
420 if (digit)
421 {
422 if(absID==nextSig)
423 {
424 // Calculate time as time of the largest digit
425 Float_t time = digit->GetTime() ;
426 Float_t aTime= digit->GetAmplitude() ;
96d5ea0d 427
a81281f7 428 // loop over input
429 Int_t nInputs = fInput;
430 if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
431 for(i = 0; i< nInputs ; i++)
432 { //loop over (possible) merge sources
433 TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
434 if(sdtclarr) {
435 Int_t sDigitEntries = sdtclarr->GetEntriesFast();
6cc75819 436
a81281f7 437 if(sDigitEntries > index[i] )
438 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
439 else
440 curSDigit = 0 ;
441
442 //May be several digits will contribute from the same input
443 while(curSDigit && (curSDigit->GetId() == absID))
444 {
445 //Shift primary to separate primaries belonging different inputs
446 Int_t primaryoffset ;
447 if(fDigInput)
448 primaryoffset = fDigInput->GetMask(i) ;
449 else
450 primaryoffset = i ;
451 curSDigit->ShiftPrimary(primaryoffset) ;
452
453 if(curSDigit->GetAmplitude()>aTime)
6cc75819 454 {
a81281f7 455 aTime = curSDigit->GetAmplitude();
456 time = curSDigit->GetTime();
6cc75819 457 }
a81281f7 458
459 *digit = *digit + *curSDigit ; //adds amplitudes of each digit
460
461 index[i]++ ;
462
463 if( sDigitEntries > index[i] )
464 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
465 else
466 curSDigit = 0 ;
467 }// while
468 }// source exists
469 }// loop over merging sources
96d5ea0d 470
96d5ea0d 471
a81281f7 472 //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
473 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
474
475 // add fluctuations for photo-electron creation
476 // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
477 Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
478 energy *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
96d5ea0d 479
a81281f7 480 //calculate and set time
481 digit->SetTime(time) ;
482
483 //Find next signal module
484 nextSig = nEMC + 1 ;
485 for(i = 0 ; i < nInputs ; i++){
486 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
487
488 if(sdigits){
489 Int_t curNext = nextSig ;
490 if(sdigits->GetEntriesFast() > index[i])
491 {
492 AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
493 if(tmpdigit)
494 {
495 curNext = tmpdigit->GetId() ;
496 }
497 }
498 if(curNext < nextSig) nextSig = curNext ;
499 }// sdigits exist
500 } // input loop
6cc75819 501
a81281f7 502 }//absID==nextSig
503
504 // add the noise now, no need for embedded events with real data
505 if(!embed)
506 energy += gRandom->Gaus(0., fPinNoise) ;
507
508
509 // JLK 26-June-2008
510 //Now digitize the energy according to the fSDigitizer method,
511 //which merely converts the energy to an integer. Later we will
512 //check that the stored value matches our allowed dynamic ranges
513 digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
514
515 //Set time resolution with final energy
516 timeResolution = GetTimeResolution(energy);
517 digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
518 AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
519 absID, energy, nextSig));
520 //Add delay to time
521 digit->SetTime(digit->GetTime()+fTimeDelay) ;
522
f21fc003 523 }// Digit pointer exists
a81281f7 524 else AliDebug(1,"Digit pointer is null");
6cc75819 525
a81281f7 526 } // for(absID = 0; absID < nEMC; absID++)
527
528 //Embed simulated amplitude (and time?) to real data digits
529 if(embed)
530 {
531 for(Int_t input = 1; input<fInput; input++)
532 {
533 TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
534 if(!realDigits)
535 {
536 AliDebug(1,"No sdigits to merge\n");
537 continue;
538 }
539
540 for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
541 {
542 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
543 if(digit2)
544 {
545 digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
546 if(digit)
547 {
548 // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
549 // and multiply to get a big int.
550 Float_t time2 = digit2->GetTime();
551 Float_t e2 = digit2->GetAmplitude();
552 CalibrateADCTime(e2,time2,digit2->GetId());
553 Float_t amp2 = fSDigitizer->Digitize(e2);
554
555 Float_t e0 = digit ->GetAmplitude();
556 if(e0>0.01)
557 AliDebug(1,Form("digit 1: Abs ID %d, amp %f, type %d, time %e; digit2: Abs ID %d, amp %f, type %d, time %e\n",
558 digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
559 digit2->GetId(),amp2, digit2->GetType(), time2 ));
560
561 // Sum amplitudes, change digit amplitude
562 digit->SetAmplitude( digit->GetAmplitude() + amp2);
563 // Rough assumption, assign time of the largest of the 2 digits,
564 // data or signal to the final digit.
565 if(amp2 > digit->GetAmplitude()) digit->SetTime(time2);
566 //Mark the new digit as embedded
567 digit->SetType(AliEMCALDigit::kEmbedded);
568
569 if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
570 AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
571 digit->GetId(), digit->GetAmplitude(), digit->GetType()));
572 }//digit2
573 }//digit2
574 }//loop digit 2
f21fc003 575 }//input loop
576 }//embed
a81281f7 577
f21fc003 578 //JLK is it better to call Clear() here?
579 delete sdigArray ; //We should not delete its contents
a81281f7 580
f21fc003 581 //remove digits below thresholds
582 // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
583 // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
584 // before merge in the same loop real data digits if available
585 Float_t energy = 0;
586 Float_t time = 0;
a81281f7 587 for(i = 0 ; i < nEMC ; i++)
588 {
f21fc003 589 digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
a81281f7 590 if(digit)
591 {
592 //Then get the energy in GeV units.
593 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
594 //Then digitize using the calibration constants of the ocdb
595 Float_t ampADC = energy;
596 DigitizeEnergyTime(ampADC, time, digit->GetId()) ;
597 if(ampADC < fDigitThreshold)
598 digits->RemoveAt(i) ;
599
f21fc003 600 }// digit exists
601 } // digit loop
a81281f7 602
f21fc003 603 digits->Compress() ;
a81281f7 604
f21fc003 605 Int_t ndigits = digits->GetEntriesFast() ;
a81281f7 606
f21fc003 607 //JLK 26-June-2008
608 //After we have done the summing and digitizing to create the
609 //digits, now we want to calibrate the resulting amplitude to match
610 //the dynamic range of our real data.
a81281f7 611 for (i = 0 ; i < ndigits ; i++)
612 {
f21fc003 613 digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
a81281f7 614 if(digit)
615 {
616 digit->SetIndexInList(i) ;
617 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
618 time = digit->GetTime();
619 Float_t ampADC = energy;
620 DigitizeEnergyTime(ampADC, time, digit->GetId());
621 digit->SetAmplitude(ampADC) ;
622 digit->SetTime(time);
623 // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
f21fc003 624 }// digit exists
625 }//Digit loop
a81281f7 626
7e1d9a9b 627 }
628 else AliFatal("EMCALLoader is NULL!") ;
6cc75819 629
814ad4bf 630}
61e0abb5 631
6cc75819 632//_____________________________________________________________________
633void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
a81281f7 634{
ab6a174f 635 // JLK 26-June-2008
e4dcc484 636 // Returns digitized value of the energy in a cell absId
ab6a174f 637 // using the calibration constants stored in the OCDB
638 // or default values if no CalibData object is found.
639 // This effectively converts everything to match the dynamic range
640 // of the real data we will collect
641 //
53f1c378 642 // Load Geometry
643 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
d62a891a 644
a81281f7 645 if (geom==0)
646 {
96d5ea0d 647 AliFatal("Did not get geometry from EMCALLoader");
6cc75819 648 return;
e4dcc484 649 }
5eb74a51 650
6cc75819 651 Int_t iSupMod = -1;
652 Int_t nModule = -1;
653 Int_t nIphi = -1;
654 Int_t nIeta = -1;
655 Int_t iphi = -1;
656 Int_t ieta = -1;
657 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
658 if(!bCell)
a81281f7 659 Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
6cc75819 660 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
96d5ea0d 661
a81281f7 662 if(fCalibData)
663 {
6cc75819 664 fADCpedestalEC = fCalibData->GetADCpedestal (iSupMod,ieta,iphi);
665 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
666 fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
031c3928 667 fTimeChannel = fCalibData->GetTimeChannel (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
a81281f7 668 fTimeChannelDecal = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
6cc75819 669 }
670
671 //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
672 energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal ;
673 time += fTimeChannel-fTimeChannelDecal;
674
a81281f7 675 if(energy > fNADCEC )
676 energy = fNADCEC ;
6cc75819 677}
829ba234 678
6cc75819 679//_____________________________________________________________________
52e22b66 680void AliEMCALDigitizer::Decalibrate(AliEMCALDigit *digit)
a81281f7 681{
682 // Decalibrate, used in Trigger digits
683 // Load Geometry
684 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
63c22917 685
a81281f7 686 if (!geom)
687 {
688 AliFatal("Did not get geometry from EMCALLoader");
689 return;
690 }
63c22917 691
a81281f7 692 Int_t absId = digit->GetId();
693 Int_t iSupMod = -1;
694 Int_t nModule = -1;
695 Int_t nIphi = -1;
696 Int_t nIeta = -1;
697 Int_t iphi = -1;
698 Int_t ieta = -1;
52e22b66 699
a81281f7 700 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
52e22b66 701
a81281f7 702 if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
703 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
63c22917 704
a81281f7 705 if (fCalibData)
706 {
707 fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
708 Float_t factor = 1./(fADCchannelEC/0.0162);
52e22b66 709
a81281f7 710 *digit = *digit * factor;
711 }
63c22917 712}
713
714//_____________________________________________________________________
6cc75819 715void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
a81281f7 716{
6cc75819 717 // Returns the energy in a cell absId with a given adc value
718 // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
719 // Used in case of embedding, transform ADC counts from real event into energy
720 // so that we can add the energy of the simulated sdigits which are in energy
721 // units.
722 // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
723 // or with time out of window
724
725 // Load Geometry
726 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
727
a81281f7 728 if (!geom)
729 {
6cc75819 730 AliFatal("Did not get geometry from EMCALLoader");
731 return;
732 }
733
734 Int_t iSupMod = -1;
735 Int_t nModule = -1;
736 Int_t nIphi = -1;
737 Int_t nIeta = -1;
738 Int_t iphi = -1;
739 Int_t ieta = -1;
740 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
a81281f7 741 if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
6cc75819 742 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
743
a81281f7 744 if(fCalibData)
745 {
6cc75819 746 fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
747 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
031c3928 748 fTimeChannel = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
6cc75819 749 }
750
751 adc = adc * fADCchannelEC - fADCpedestalEC;
752 time -= fTimeChannel;
96d5ea0d 753
61e0abb5 754}
755
6cc75819 756
61e0abb5 757//____________________________________________________________________________
f21fc003 758void AliEMCALDigitizer::Digitize(Option_t *option)
05a92d59 759{
45fa49ca 760 // Steering method to process digitization for events
761 // in the range from fFirstEvent to fLastEvent.
762 // This range is optionally set by SetEventRange().
763 // if fLastEvent=-1, then process events until the end.
764 // by default fLastEvent = fFirstEvent (process only one event)
05a92d59 765
88cb7938 766 if (!fInit) { // to prevent overwrite existing file
f21fc003 767 Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
88cb7938 768 return ;
769 }
770
814ad4bf 771 if (strstr(option,"print")) {
b42d4197 772
88cb7938 773 Print();
814ad4bf 774 return ;
775 }
776
61e0abb5 777 if(strstr(option,"tim"))
778 gBenchmark->Start("EMCALDigitizer");
5dee926e 779
33c3c91a 780 AliRunLoader *rl = AliRunLoader::Instance();
5dee926e 781 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
96d5ea0d 782 Int_t nEvents = 0;
783 if(!emcalLoader){
784 AliFatal("Did not get the Loader");
30bc5594 785 }
96d5ea0d 786 else{
64942713 787
96d5ea0d 788 if (fLastEvent == -1) {
789 fLastEvent = rl->GetNumberOfEvents() - 1 ;
790 }
f21fc003 791 else if (fDigInput)
96d5ea0d 792 fLastEvent = fFirstEvent ; // what is this ??
793
794 nEvents = fLastEvent - fFirstEvent + 1;
795 Int_t ievent = -1;
88cb7938 796
8cc543cb 797 AliEMCALGeometry *geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
798 const Int_t nTRU = geom->GetNTotalTRU();
799 TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit", nTRU*96);
800 TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
88cb7938 801
96d5ea0d 802 rl->LoadSDigits("EMCAL");
803 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
804
805 rl->GetEvent(ievent);
806
807 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
808
809 WriteDigits() ;
810
811 //Trigger Digits
812 //-------------------------------------
813
88cb7938 814
96d5ea0d 815 Digits2FastOR(digitsTMP, digitsTRG);
816
817 WriteDigits(digitsTRG);
818
819 (emcalLoader->TreeD())->Fill();
820
821 emcalLoader->WriteDigits( "OVERWRITE");
96d5ea0d 822
823 Unload();
824
825 digitsTRG ->Delete();
826 digitsTMP ->Delete();
827
828 //-------------------------------------
829
830 if(strstr(option,"deb"))
831 PrintDigits(option);
832 if(strstr(option,"table")) gObjectTable->Print();
833
834 //increment the total number of Digits per run
835 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
836 }//loop
f21fc003 837
96d5ea0d 838 }//loader exists
64942713 839
61e0abb5 840 if(strstr(option,"tim")){
841 gBenchmark->Stop("EMCALDigitizer");
0bd264d5 842 Float_t cputime = gBenchmark->GetCpuTime("EMCALDigitizer");
843 Float_t avcputime = cputime;
844 if(nEvents==0) avcputime = 0 ;
845 AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
88cb7938 846 }
61e0abb5 847}
848
a2f8e711 849//__________________________________________________________________
850Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
851{
a81281f7 852 // Assign a smeared time to the digit, from observed distribution in LED system (?).
a2f8e711 853 // From F. Blanco
854 Float_t res = -1;
a81281f7 855 if (energy > 0)
856 {
a2f8e711 857 res = TMath::Sqrt(fTimeResolutionPar0 +
858 fTimeResolutionPar1/(energy*energy) );
859 }
b2a891e0 860 // parametrization above is for ns. Convert to seconds:
861 return res*1e-9;
a2f8e711 862}
863
a81281f7 864//____________________________________________________________________________
916f1e76 865void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
866{
a81281f7 867 // FEE digits afterburner to produce TRG digits
de39a0ff 868 // we are only interested in the FEE digit deposited energy
869 // to be converted later into a voltage value
870
871 // push the FEE digit to its associated FastOR (numbered from 0:95)
872 // TRU is in charge of summing module digits
873
874 AliRunLoader *runLoader = AliRunLoader::Instance();
875
de39a0ff 876 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
6995e8b7 877 if (!emcalLoader) AliFatal("Did not get the Loader");
878
a81281f7 879 const AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
880
881 // build FOR from simulated digits
882 // and xfer to the corresponding TRU input (mapping)
883
884 TClonesArray* sdigits = emcalLoader->SDigits();
885
886 TClonesArray *digits = (TClonesArray*)sdigits->Clone();
887
888 AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
889
890 TIter NextDigit(digits);
891 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
892 {
893 if (IsDead(digit)) continue;
894
895 Decalibrate(digit);
896
897 Int_t id = digit->GetId();
898
899 Int_t trgid;
900 if (geom && geom->GetFastORIndexFromCellIndex(id , trgid))
901 {
902 AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
7e1d9a9b 903
a81281f7 904 AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
7d525d9d 905
a81281f7 906 if (!d)
96d5ea0d 907 {
a81281f7 908 new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
909 d = (AliEMCALDigit*)digitsTMP->At(trgid);
910 d->SetId(trgid);
911 }
912 else
913 {
914 *d = *d + *digit;
96d5ea0d 915 }
a81281f7 916 }
917 }
918
919 if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
920
921 Int_t nSamples = geom->GetNTotalTRU();
922 Int_t *timeSamples = new Int_t[nSamples];
923
924 NextDigit = TIter(digitsTMP);
925 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
926 {
927 if (digit)
928 {
929 Int_t id = digit->GetId();
930 Float_t time = 50.e-9;
7e1d9a9b 931
a81281f7 932 Double_t depositedEnergy = 0.;
933 for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
7e1d9a9b 934
a81281f7 935 if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
936
937 // FIXME: Check digit time!
938 if (depositedEnergy) {
939 depositedEnergy += gRandom->Gaus(0., .08);
940 DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
941
942 for (Int_t j=0;j<nSamples;j++) {
943 if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
944 timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
7e1d9a9b 945 }
a81281f7 946
947 new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
948
949 if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
96d5ea0d 950 }
a81281f7 951 }
952 }
953
954 delete [] timeSamples;
955
956 if (digits && digits->GetEntriesFast()) digits->Delete();
916f1e76 957}
958
a81281f7 959//____________________________________________________________________________
916f1e76 960void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
961{
a81281f7 962 // parameters:
963 // id: 0..95
964 const Int_t reso = 12; // 11-bit resolution ADC
965 const Double_t vFSR = 2.; // Full scale input voltage range 2V (p-p)
966//const Double_t dNe = 125; // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
967 const Double_t dNe = 125/1.3; // F-ALTRO max V. FEE: factor 4
968 const Double_t vA = .136e-6; // CSP output range: 0.136uV/e-
969 const Double_t rise = 50e-9; // rise time (10-90%) of the FastOR signal before shaping
916f1e76 970
a81281f7 971 const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
916f1e76 972
a81281f7 973 Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
974
975 TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
976 signalF.SetParameter( 0, vV );
977 signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
978 signalF.SetParameter( 2, rise );
916f1e76 979
a81281f7 980 for (Int_t iTime=0; iTime<nSamples; iTime++)
981 {
982 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
983 // probably plan an access to OCDB
28bcaaaa 984 Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
985 if (TMath::Abs(sig) > vFSR/2.) {
986 AliError("Signal overflow!");
987 timeSamples[iTime] = (1 << reso) - 1;
988 } else {
989 AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
990 timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
991 }
a81281f7 992 }
916f1e76 993}
994
05a92d59 995//____________________________________________________________________________
996Bool_t AliEMCALDigitizer::Init()
997{
998 // Makes all memory allocations
88cb7938 999 fInit = kTRUE ;
33c3c91a 1000 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
61d80ce0 1001
5dee926e 1002 if ( emcalLoader == 0 ) {
1003 Fatal("Init", "Could not obtain the AliEMCALLoader");
05a92d59 1004 return kFALSE;
1005 }
61d80ce0 1006
45fa49ca 1007 fFirstEvent = 0 ;
1008 fLastEvent = fFirstEvent ;
1009
f21fc003 1010 if(fDigInput)
1011 fInput = fDigInput->GetNinputs() ;
88cb7938 1012 else
1013 fInput = 1 ;
1014
1015 fInputFileNames = new TString[fInput] ;
1016 fEventNames = new TString[fInput] ;
1017 fInputFileNames[0] = GetTitle() ;
1018 fEventNames[0] = fEventFolderName.Data() ;
1019 Int_t index ;
1020 for (index = 1 ; index < fInput ; index++) {
f21fc003 1021 fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1022 TString tempo = fDigInput->GetInputFolderName(index) ;
1023 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
05a92d59 1024 }
f21fc003 1025
53f1c378 1026 //Calibration instance
1027 fCalibData = emcalLoader->CalibData();
88cb7938 1028 return fInit ;
05a92d59 1029}
1030
1031//____________________________________________________________________________
1032void AliEMCALDigitizer::InitParameters()
14ce0a6e 1033{
29b5d492 1034 // Parameter initialization for digitizer
6569f329 1035
1036 // Get the parameters from the OCDB via the loader
1037 AliRunLoader *rl = AliRunLoader::Instance();
1038 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1039 AliEMCALSimParam * simParam = 0x0;
1040 if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1041
1042 if(!simParam){
61d80ce0 1043 simParam = AliEMCALSimParam::GetInstance();
1044 AliWarning("Simulation Parameters not available in OCDB?");
6569f329 1045 }
61d80ce0 1046
5fb45456 1047 fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400; // electrons per GeV
1048 fGainFluctuations = simParam->GetGainFluctuations() ; // 15.0;
1049
6569f329 1050 fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
33a52e55 1051 if (fPinNoise < 0.0001 )
88cb7938 1052 Warning("InitParameters", "No noise added\n") ;
6cc75819 1053 fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
6569f329 1054 fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
a2f8e711 1055 fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1056 fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
f0a6dc6f 1057 fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
61d80ce0 1058
97075cd8 1059 // These defaults are normally not used.
1060 // Values are read from calibration database instead
33a52e55 1061 fADCchannelEC = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
6cc75819 1062 fADCpedestalEC = 0.0 ; // GeV
1063 fADCchannelECDecal = 1.0; // No decalibration by default
1064 fTimeChannel = 0.0; // No time calibration by default
1065 fTimeChannelDecal = 0.0; // No time decalibration by default
33a52e55 1066
6cc75819 1067 fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
61d80ce0 1068
5fb45456 1069 AliDebug(2,Form("Mean Photon Electron %d, Gain Fluct. %2.1f; Noise: APD %f, Time %f; Digit Threshold %d,Time Resolution Par0 %g Par1 %g,NADCEC %d",
1070 fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1071
05a92d59 1072}
61e0abb5 1073
61e0abb5 1074//__________________________________________________________________
1963b290 1075void AliEMCALDigitizer::Print1(Option_t * option)
6cc75819 1076{ // 19-nov-04 - just for convenience
1963b290 1077 Print();
1078 PrintDigits(option);
1079}
1080
61e0abb5 1081//__________________________________________________________________
6cc75819 1082void AliEMCALDigitizer::Print (Option_t * ) const
88cb7938 1083{
1084 // Print Digitizer's parameters
fdebddeb 1085 printf("Print: \n------------------- %s -------------", GetName() ) ;
88cb7938 1086 if( strcmp(fEventFolderName.Data(), "") != 0 ){
1087 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1088
1089 Int_t nStreams ;
f21fc003 1090 if (fDigInput)
88cb7938 1091 nStreams = GetNInputStreams() ;
1092 else
1093 nStreams = fInput ;
61d80ce0 1094
53e430a3 1095 AliRunLoader *rl=0;
88cb7938 1096
1097 Int_t index = 0 ;
1098 for (index = 0 ; index < nStreams ; index++) {
1099 TString tempo(fEventNames[index]) ;
1100 tempo += index ;
5dee926e 1101 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
96d5ea0d 1102 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
5dee926e 1103 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
96d5ea0d 1104 if(emcalLoader){
1105 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1106 if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1107 fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
1108 printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1109 }//loader
88cb7938 1110 }
61d80ce0 1111
33c3c91a 1112 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
61d80ce0 1113
96d5ea0d 1114 if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1115 else printf("\nNULL LOADER");
61d80ce0 1116
88cb7938 1117 printf("\nWith following parameters:\n") ;
6cc75819 1118 printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
4354827d 1119 printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
88cb7938 1120 printf("---------------------------------------------------\n") ;
ea5f3389 1121 }
1122 else
fdebddeb 1123 printf("Print: AliEMCALDigitizer not initialized") ;
61e0abb5 1124}
814ad4bf 1125
61e0abb5 1126//__________________________________________________________________
14ce0a6e 1127void AliEMCALDigitizer::PrintDigits(Option_t * option)
1128{
1129 //utility method for printing digit information
61d80ce0 1130
33c3c91a 1131 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
96d5ea0d 1132 if(emcalLoader){
1133 TClonesArray * digits = emcalLoader->Digits() ;
1134 TClonesArray * sdigits = emcalLoader->SDigits() ;
1135
1136 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1137 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1138
1139 if(strstr(option,"all")){
1140 //loop over digits
1141 AliEMCALDigit * digit;
1142 printf("\nEMC digits (with primaries):\n") ;
1143 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1144 Int_t index ;
1145 for (index = 0 ; index < digits->GetEntries() ; index++) {
1146 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1147 if(digit){
1148 printf("\n%6d %8f %6.5e %4d %2d : ",
1149 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1150 Int_t iprimary;
1151 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1152 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1153 }
1154 }// digit exists
1155 }// loop
1156 }
1157 printf("\n");
1158 }// loader exists
1159 else printf("NULL LOADER, cannot print\n");
61e0abb5 1160}
05a92d59 1161
1162//__________________________________________________________________
1163Float_t AliEMCALDigitizer::TimeOfNoise(void)
04f0bda3 1164{
1165 // Calculates the time signal generated by noise
6cc75819 1166 //printf("Time noise %e\n",fTimeNoise);
1167 return gRandom->Rndm() * fTimeNoise;
05a92d59 1168}
1169
88cb7938 1170//__________________________________________________________________
1171void AliEMCALDigitizer::Unload()
1172{
fdebddeb 1173 // Unloads the SDigits and Digits
5dee926e 1174 AliRunLoader *rl=0;
1175
88cb7938 1176 Int_t i ;
1177 for(i = 1 ; i < fInput ; i++){
1178 TString tempo(fEventNames[i]) ;
5dee926e 1179 tempo += i;
1180 if ((rl = AliRunLoader::GetRunLoader(tempo)))
1181 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
88cb7938 1182 }
33c3c91a 1183 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
96d5ea0d 1184 if(emcalLoader)emcalLoader->UnloadDigits() ;
1185 else AliFatal("Did not get the loader");
88cb7938 1186}
1187
814ad4bf 1188//_________________________________________________________________________________________
9e5d2067 1189void AliEMCALDigitizer::WriteDigits()
a2f8e711 1190{ // Makes TreeD in the output file.
814ad4bf 1191 // Check if branch already exists:
1192 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1193 // already existing branches.
1194 // else creates branch with Digits, named "EMCAL", title "...",
1195 // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1196 // and names of files, from which digits are made.
96d5ea0d 1197
33c3c91a 1198 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
96d5ea0d 1199
1200 if(emcalLoader){
1201 const TClonesArray * digits = emcalLoader->Digits() ;
1202 TTree * treeD = emcalLoader->TreeD();
1203 if ( !treeD ) {
1204 emcalLoader->MakeDigitsContainer();
1205 treeD = emcalLoader->TreeD();
1206 }
1207
1208 // -- create Digits branch
1209 Int_t bufferSize = 32000 ;
1210 TBranch * digitsBranch = 0;
1211 if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1212 digitsBranch->SetAddress(&digits);
1213 AliWarning("Digits Branch already exists. Not all digits will be visible");
1214 }
1215 else
1216 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1217 //digitsBranch->SetTitle(fEventFolderName);
1218
1219 // treeD->Fill() ;
1220 /*
1221 emcalLoader->WriteDigits("OVERWRITE");
1222 emcalLoader->WriteDigitizer("OVERWRITE");
1223
1224 Unload() ;
1225 */
1226
1227 }// loader exists
1228 else AliFatal("Loader not available");
916f1e76 1229}
88cb7938 1230
916f1e76 1231//__________________________________________________________________
1232void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
a2f8e711 1233{ // overloaded method
61d80ce0 1234 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1235 if(emcalLoader){
96d5ea0d 1236
1237 TTree* treeD = emcalLoader->TreeD();
1238 if (!treeD)
61d80ce0 1239 {
1240 emcalLoader->MakeDigitsContainer();
1241 treeD = emcalLoader->TreeD();
1242 }
96d5ea0d 1243
1244 // -- create Digits branch
1245 Int_t bufferSize = 32000;
1246
1247 if (TBranch* triggerBranch = treeD->GetBranch(branchName))
61d80ce0 1248 {
1249 triggerBranch->SetAddress(&digits);
1250 }
96d5ea0d 1251 else
61d80ce0 1252 {
1253 treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1254 }
96d5ea0d 1255
1256 // treeD->Fill();
1257 }// loader exists
1258 else AliFatal("Loader not available");
814ad4bf 1259}
61e0abb5 1260
6995e8b7 1261//__________________________________________________________________
1262Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1263{
1264 AliRunLoader *runLoader = AliRunLoader::Instance();
1265 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1266 if (!emcalLoader) AliFatal("Did not get the Loader");
1267
1268 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1269 if (!caloPed) {
1270 AliWarning("Could not access pedestal data! No dead channel removal applied");
1271 return kFALSE;
1272 }
1273
1274 // Load Geometry
1275 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1276 if (!geom) AliFatal("Did not get geometry from EMCALLoader");
1277
1278 Int_t absId = digit->GetId();
1279 Int_t iSupMod = -1;
1280 Int_t nModule = -1;
1281 Int_t nIphi = -1;
1282 Int_t nIeta = -1;
1283 Int_t iphi = -1;
1284 Int_t ieta = -1;
1285
1286 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1287
1288 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1289 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1290
1291 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1292
1293 if (channelStatus == AliCaloCalibPedestal::kDead)
1294 return kTRUE;
1295 else
1296 return kFALSE;
1297}
1298
1299
1300//__________________________________________________________________
1301Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1302{
1303 AliRunLoader *runLoader = AliRunLoader::Instance();
1304 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1305 if (!emcalLoader) AliFatal("Did not get the Loader");
1306
1307 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1308 if (!caloPed) {
1309 AliWarning("Could not access pedestal data! No dead channel removal applied");
1310 return kFALSE;
1311 }
1312
1313 // Load Geometry
1314 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1315 if (!geom) AliFatal("Did not get geometry from EMCALLoader");
1316
1317 Int_t iSupMod = -1;
1318 Int_t nModule = -1;
1319 Int_t nIphi = -1;
1320 Int_t nIeta = -1;
1321 Int_t iphi = -1;
1322 Int_t ieta = -1;
1323
1324 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1325
1326 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1327 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1328
1329 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1330
1331 if (channelStatus == AliCaloCalibPedestal::kDead)
1332 return kTRUE;
1333 else
1334 return kFALSE;
1335}
1336