1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //_________________________________________________________________________
20 //////////////////////////////////////////////////////////////////////////////
21 // Class performs digitization of Summable digits from simulated data
23 // In addition it performs mixing/embedding of summable digits from different events.
25 // For each event 3 branches are created in TreeD:
26 // "EMCAL" - list of digits
27 // "EMCALTRG" - list of trigger digits
28 // "AliEMCALDigitizer" - AliEMCALDigitizer with all parameters used in digitization
31 ////////////////////////////////////////////////////////////////////////////////////
33 //*-- Author: Sahal Yacoob (LBL)
34 // based on : AliEMCALDigitizer
36 // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
37 // of new IO (a la PHOS)
38 // November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry
39 // July 2011 GCB: Digitizer modified to accomodate embedding.
40 // Time calibration added. Decalibration possibility of energy and time added
41 //_________________________________________________________________________________
43 // --- ROOT system ---
47 #include <TBenchmark.h>
49 #include <TObjectTable.h>
54 // --- AliRoot header files ---
57 #include "AliDigitizationInput.h"
58 #include "AliRunLoader.h"
59 #include "AliCDBManager.h"
60 #include "AliCDBEntry.h"
61 #include "AliEMCALDigit.h"
63 #include "AliEMCALLoader.h"
64 #include "AliEMCALDigitizer.h"
65 #include "AliEMCALSDigitizer.h"
66 #include "AliEMCALGeometry.h"
67 #include "AliEMCALCalibData.h"
68 #include "AliEMCALSimParam.h"
69 #include "AliEMCALRawDigit.h"
73 Double_t HeavisideTheta(Double_t x)
77 if (x > 0.) signal = 1.;
82 Double_t AnalogFastORFunction(Double_t *x, Double_t *par)
95 return (((0.8*(-((TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-t + t0)/(C1*R1))*
96 TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2)) +
97 C1*C2*R1*R2*(1 - (C2*TMath::Power(TMath::E(),(-t + t0)/(C2*R2))*R2)/(-(C1*R1) + C2*R2)))*v0*
98 HeavisideTheta(t - t0))/tr
100 (TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C1*R1))*
101 TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2) +
102 (C1*TMath::Power(C2,2)*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C2*R2))*
103 R1*TMath::Power(R2,2))/(C1*R1 - C2*R2))*v0*
104 HeavisideTheta(t - t0 - 1.25*tr))/tr)/(C2*R1));
108 ClassImp(AliEMCALDigitizer)
111 //____________________________________________________________________________
112 AliEMCALDigitizer::AliEMCALDigitizer()
113 : AliDigitizer("",""),
118 fInputFileNames(0x0),
121 fMeanPhotonElectron(0),
122 fGainFluctuations(0),
126 fTimeResolutionPar0(0),
127 fTimeResolutionPar1(0),
130 fADCchannelECDecal(0),
132 fTimeChannelDecal(0),
134 fEventFolderName(""),
142 fDigInput = 0 ; // We work in the standalone mode
145 //____________________________________________________________________________
146 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
147 : AliDigitizer("EMCALDigitizer", alirunFileName),
148 fDefaultInit(kFALSE),
155 fMeanPhotonElectron(0),
156 fGainFluctuations(0),
160 fTimeResolutionPar0(0),
161 fTimeResolutionPar1(0),
164 fADCchannelECDecal(0),
166 fTimeChannelDecal(0),
168 fEventFolderName(eventFolderName),
177 fDigInput = 0 ; // We work in the standalone mode
180 //____________________________________________________________________________
181 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
182 : AliDigitizer(d.GetName(),d.GetTitle()),
183 fDefaultInit(d.fDefaultInit),
184 fDigitsInRun(d.fDigitsInRun),
187 fInputFileNames(d.fInputFileNames),
188 fEventNames(d.fEventNames),
189 fDigitThreshold(d.fDigitThreshold),
190 fMeanPhotonElectron(d.fMeanPhotonElectron),
191 fGainFluctuations(d.fGainFluctuations),
192 fPinNoise(d.fPinNoise),
193 fTimeNoise(d.fTimeNoise),
194 fTimeDelay(d.fTimeDelay),
195 fTimeResolutionPar0(d.fTimeResolutionPar0),
196 fTimeResolutionPar1(d.fTimeResolutionPar1),
197 fADCchannelEC(d.fADCchannelEC),
198 fADCpedestalEC(d.fADCpedestalEC),
199 fADCchannelECDecal(d.fADCchannelECDecal),
200 fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
202 fEventFolderName(d.fEventFolderName),
203 fFirstEvent(d.fFirstEvent),
204 fLastEvent(d.fLastEvent),
205 fCalibData(d.fCalibData),
206 fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
211 //____________________________________________________________________________
212 AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
213 : AliDigitizer(rd,"EMCALDigitizer"),
214 fDefaultInit(kFALSE),
221 fMeanPhotonElectron(0),
222 fGainFluctuations(0),
226 fTimeResolutionPar0(0.),
227 fTimeResolutionPar1(0.),
230 fADCchannelECDecal(0),
232 fTimeChannelDecal(0),
240 // ctor Init() is called by RunDigitizer
242 fEventFolderName = fDigInput->GetInputFolderName(0) ;
243 SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
247 //____________________________________________________________________________
248 AliEMCALDigitizer::~AliEMCALDigitizer()
251 delete [] fInputFileNames ;
252 delete [] fEventNames ;
253 if (fSDigitizer) delete fSDigitizer;
256 //____________________________________________________________________________
257 void AliEMCALDigitizer::Digitize(Int_t event)
260 // Makes the digitization of the collected summable digits
261 // for this it first creates the array of all EMCAL modules
262 // filled with noise and after that adds contributions from
263 // SDigits. This design helps to avoid scanning over the
264 // list of digits to add contribution of any new SDigit.
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
271 static int nEMC=0; //max number of digits possible
272 AliRunLoader *rl = AliRunLoader::Instance();
273 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
276 Int_t readEvent = event ;
277 // fDigInput is data member from AliDigitizer
279 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
280 AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
281 readEvent, GetTitle(), fEventFolderName.Data())) ;
282 rl->GetEvent(readEvent);
284 TClonesArray * digits = emcalLoader->Digits() ;
285 digits->Delete() ; //correct way to clear array when memory is
286 //allocated by objects stored in array
289 AliEMCALGeometry *geom = 0;
290 if (!rl->GetAliRun()) {
291 AliFatal("Could not get AliRun from runLoader");
294 AliEMCAL * emcal = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
295 geom = emcal->GetGeometry();
296 nEMC = geom->GetNCells();
297 AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
302 digits->Expand(nEMC) ;
304 // RS create a digitizer on the fly
305 if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
306 fSDigitizer->SetEventRange(0, -1) ;
308 //take all the inputs to add together and load the SDigits
309 TObjArray * sdigArray = new TObjArray(fInput) ;
310 sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
312 Int_t embed = kFALSE;
313 for(i = 1 ; i < fInput ; i++){
314 TString tempo(fEventNames[i]) ;
317 AliRunLoader * rl2 = AliRunLoader::GetRunLoader(tempo) ;
320 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
323 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
324 Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
326 // rl2->LoadDigits();
327 rl2->GetEvent(readEvent);
328 if(rl2->GetDetectorLoader("EMCAL"))
330 AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
333 //Merge 2 simulated sdigits
334 if(emcalLoader2->SDigits()){
335 TClonesArray* sdigits2 = emcalLoader2->SDigits();
336 sdigArray->AddAt(sdigits2, i) ;
337 // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
338 // do not smear energy of embedded, do not add noise to any sdigits
339 if(sdigits2->GetEntriesFast()>0){
340 //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
341 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
342 if(digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded){
349 else AliFatal("EMCAL Loader of second event not available!");
352 else Fatal("Digitize", "Loader of second input not found");
355 //Find the first tower with signal
356 Int_t nextSig = nEMC + 1 ;
357 TClonesArray * sdigits ;
358 for(i = 0 ; i < fInput ; i++){
359 if(i > 0 && embed) continue;
360 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
362 if (sdigits->GetEntriesFast() ){
363 AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
365 Int_t curNext = sd->GetId() ;
366 if(curNext < nextSig)
368 AliDebug(1,Form("input %i : #sdigits %i \n",
369 i, sdigits->GetEntriesFast()));
371 }//First entry is not NULL
373 Info("Digitize", "NULL sdigit pointer");
376 }//There is at least one entry
378 Info("Digitize", "NULL sdigits array");
381 }// SDigits array exists
383 Info("Digitizer","Null sdigit pointer");
387 AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
389 TArrayI index(fInput) ;
390 index.Reset() ; //Set all indexes to zero
392 AliEMCALDigit * digit ;
393 AliEMCALDigit * curSDigit ;
395 //Put Noise contribution, smear time and energy
396 Float_t timeResolution = 0;
397 for(absID = 0; absID < nEMC; absID++){ // Nov 30, 2006 by PAI; was from 1 to nEMC
400 // amplitude set to zero, noise will be added later
401 Float_t noiseTime = 0.;
402 if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
403 new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
404 //look if we have to add signal?
405 digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
411 // Calculate time as time of the largest digit
412 Float_t time = digit->GetTime() ;
413 Float_t aTime= digit->GetAmplitude() ;
416 Int_t nInputs = fInput;
417 if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
418 for(i = 0; i< nInputs ; i++){ //loop over (possible) merge sources
419 TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
421 Int_t sDigitEntries = sdtclarr->GetEntriesFast();
423 if(sDigitEntries > index[i] )
424 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
428 //May be several digits will contribute from the same input
429 while(curSDigit && (curSDigit->GetId() == absID)){
430 //Shift primary to separate primaries belonging different inputs
431 Int_t primaryoffset ;
433 primaryoffset = fDigInput->GetMask(i) ;
436 curSDigit->ShiftPrimary(primaryoffset) ;
438 if(curSDigit->GetAmplitude()>aTime) {
439 aTime = curSDigit->GetAmplitude();
440 time = curSDigit->GetTime();
443 *digit = *digit + *curSDigit ; //adds amplitudes of each digit
447 if( sDigitEntries > index[i] )
448 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
453 }// loop over merging sources
456 //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
457 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
459 // add fluctuations for photo-electron creation
460 // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
461 Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
462 energy *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
464 //calculate and set time
465 digit->SetTime(time) ;
467 //Find next signal module
469 for(i = 0 ; i < nInputs ; i++){
470 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
473 Int_t curNext = nextSig ;
474 if(sdigits->GetEntriesFast() > index[i])
476 AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
479 curNext = tmpdigit->GetId() ;
482 if(curNext < nextSig) nextSig = curNext ;
488 // add the noise now, no need for embedded events with real data
490 energy += gRandom->Gaus(0., fPinNoise) ;
494 //Now digitize the energy according to the fSDigitizer method,
495 //which merely converts the energy to an integer. Later we will
496 //check that the stored value matches our allowed dynamic ranges
497 digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
499 //Set time resolution with final energy
500 timeResolution = GetTimeResolution(energy);
501 digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
502 AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
503 absID, energy, nextSig));
505 digit->SetTime(digit->GetTime()+fTimeDelay) ;
507 }// Digit pointer exists
509 Info("Digitizer","Digit pointer is null");
511 } // for(absID = 0; absID < nEMC; absID++)
516 //Embed simulated amplitude (and time?) to real data digits
518 for(Int_t input = 1; input<fInput; input++){
519 TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
521 AliDebug(1,"No sdigits to merge\n");
524 for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++){
525 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
527 digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
530 // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
531 // and multiply to get a big int.
532 Float_t time2 = digit2->GetTime();
533 Float_t e2 = digit2->GetAmplitude();
534 CalibrateADCTime(e2,time2,digit2->GetId());
535 Float_t amp2 = fSDigitizer->Digitize(e2);
537 Float_t e0 = digit ->GetAmplitude();
539 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",
540 digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
541 digit2->GetId(),amp2, digit2->GetType(), time2 ));
543 // Sum amplitudes, change digit amplitude
544 digit->SetAmplitude( digit->GetAmplitude() + amp2);
545 // Rough assumption, assign time of the largest of the 2 digits,
546 // data or signal to the final digit.
547 if(amp2 > digit->GetAmplitude()) digit->SetTime(time2);
548 //Mark the new digit as embedded
549 digit->SetType(AliEMCALDigit::kEmbedded);
551 if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
552 AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
553 digit->GetId(), digit->GetAmplitude(), digit->GetType()));
560 //JLK is it better to call Clear() here?
561 delete sdigArray ; //We should not delete its contents
563 //remove digits below thresholds
564 // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
565 // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
566 // before merge in the same loop real data digits if available
569 for(i = 0 ; i < nEMC ; i++){
570 digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
573 //Then get the energy in GeV units.
574 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
575 //Then digitize using the calibration constants of the ocdb
576 Float_t ampADC = energy;
577 DigitizeEnergyTime(ampADC, time, digit->GetId()) ;
578 if(ampADC < fDigitThreshold)
579 digits->RemoveAt(i) ;
586 Int_t ndigits = digits->GetEntriesFast() ;
589 //After we have done the summing and digitizing to create the
590 //digits, now we want to calibrate the resulting amplitude to match
591 //the dynamic range of our real data.
592 for (i = 0 ; i < ndigits ; i++) {
593 digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
595 digit->SetIndexInList(i) ;
596 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
597 time = digit->GetTime();
598 Float_t ampADC = energy;
599 DigitizeEnergyTime(ampADC, time, digit->GetId());
600 digit->SetAmplitude(ampADC) ;
601 digit->SetTime(time);
602 // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
607 else AliFatal("EMCALLoader is NULL!") ;
611 //_____________________________________________________________________
612 void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
615 // Returns digitized value of the energy in a cell absId
616 // using the calibration constants stored in the OCDB
617 // or default values if no CalibData object is found.
618 // This effectively converts everything to match the dynamic range
619 // of the real data we will collect
622 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
625 AliFatal("Did not get geometry from EMCALLoader");
635 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
637 Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
638 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
641 fADCpedestalEC = fCalibData->GetADCpedestal (iSupMod,ieta,iphi);
642 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
643 fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
644 fTimeChannel = fCalibData->GetTimeChannel (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
645 fTimeChannelDecal = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
648 //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
649 energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal ;
650 time += fTimeChannel-fTimeChannelDecal;
652 if(energy > fNADCEC )
656 //_____________________________________________________________________
657 void AliEMCALDigitizer::DecalibrateEnergy(Double_t& energy, const Int_t absId)
660 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
663 AliFatal("Did not get geometry from EMCALLoader");
673 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
675 Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
676 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
679 fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
682 energy /= fADCchannelEC/0.0162;
685 //_____________________________________________________________________
686 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
688 // Returns the energy in a cell absId with a given adc value
689 // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
690 // Used in case of embedding, transform ADC counts from real event into energy
691 // so that we can add the energy of the simulated sdigits which are in energy
693 // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
694 // or with time out of window
697 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
700 AliFatal("Did not get geometry from EMCALLoader");
710 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
712 Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
713 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
716 fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
717 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
718 fTimeChannel = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
721 adc = adc * fADCchannelEC - fADCpedestalEC;
722 time -= fTimeChannel;
728 //____________________________________________________________________________
729 void AliEMCALDigitizer::Digitize(Option_t *option)
731 // Steering method to process digitization for events
732 // in the range from fFirstEvent to fLastEvent.
733 // This range is optionally set by SetEventRange().
734 // if fLastEvent=-1, then process events until the end.
735 // by default fLastEvent = fFirstEvent (process only one event)
737 if (!fInit) { // to prevent overwrite existing file
738 Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
742 if (strstr(option,"print")) {
748 if(strstr(option,"tim"))
749 gBenchmark->Start("EMCALDigitizer");
751 AliRunLoader *rl = AliRunLoader::Instance();
752 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
755 AliFatal("Did not get the Loader");
759 if (fLastEvent == -1) {
760 fLastEvent = rl->GetNumberOfEvents() - 1 ;
763 fLastEvent = fFirstEvent ; // what is this ??
765 nEvents = fLastEvent - fFirstEvent + 1;
768 TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit", 32*96);
769 TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", 32*96);
771 rl->LoadSDigits("EMCAL");
772 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
774 rl->GetEvent(ievent);
776 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
781 //-------------------------------------
784 Digits2FastOR(digitsTMP, digitsTRG);
786 WriteDigits(digitsTRG);
788 (emcalLoader->TreeD())->Fill();
790 emcalLoader->WriteDigits( "OVERWRITE");
794 digitsTRG ->Delete();
795 digitsTMP ->Delete();
797 //-------------------------------------
799 if(strstr(option,"deb"))
801 if(strstr(option,"table")) gObjectTable->Print();
803 //increment the total number of Digits per run
804 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
809 if(strstr(option,"tim")){
810 gBenchmark->Stop("EMCALDigitizer");
811 Float_t cputime = gBenchmark->GetCpuTime("EMCALDigitizer");
812 Float_t avcputime = cputime;
813 if(nEvents==0) avcputime = 0 ;
814 AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
818 //__________________________________________________________________
819 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
824 res = TMath::Sqrt(fTimeResolutionPar0 +
825 fTimeResolutionPar1/(energy*energy) );
827 // parametrization above is for ns. Convert to seconds:
831 //____________________________________________________________________________
832 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
834 // FEE digits afterburner to produce TRG digits
835 // we are only interested in the FEE digit deposited energy
836 // to be converted later into a voltage value
838 // push the FEE digit to its associated FastOR (numbered from 0:95)
839 // TRU is in charge of summing module digits
841 AliRunLoader *runLoader = AliRunLoader::Instance();
843 AliRun* run = runLoader->GetAliRun();
845 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
847 AliFatal("Did not get the Loader");
850 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(run->GetDetector("EMCAL"));
852 AliEMCALGeometry* geom = emcal->GetGeometry();
854 // build FOR from simulated digits
855 // and xfer to the corresponding TRU input (mapping)
857 TClonesArray* digits = emcalLoader->SDigits();
859 AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
861 TIter NextDigit(digits);
862 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
864 Int_t id = digit->GetId();
867 if (geom && geom->GetFastORIndexFromCellIndex(id , trgid))
869 AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
871 AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
875 new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
876 d = (AliEMCALDigit*)digitsTMP->At(trgid);
886 if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
888 TF1 g4noiseF("g4noise","[0]*exp(-0.5*((x-[1])/[2])**4)",-10,10);
889 g4noiseF.SetParameter( 0, 1 ); //
890 g4noiseF.SetParameter( 1, 0 ); //
891 g4noiseF.SetParameter( 2, .2 ); //
894 Int_t *timeSamples = new Int_t[nSamples];
896 NextDigit = TIter(digitsTMP);
897 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
901 Int_t id = digit->GetId();
902 Float_t time = 50.e-9;
904 Double_t depositedEnergy = 0.;
905 for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
907 if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
909 // FIXME: Check digit time!
912 // depositedEnergy += depositedEnergy * g4noiseF.GetRandom();
914 DecalibrateEnergy(depositedEnergy, id);
916 DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
918 for (Int_t j=0;j<nSamples;j++)
920 timeSamples[j] = ((j << 12) & 0xFF000) | (timeSamples[j] & 0xFFF);
923 new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
925 if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
930 delete [] timeSamples;
932 else AliFatal("Could not get AliEMCAL");
937 //____________________________________________________________________________
938 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
942 const Int_t reso = 11; // 11-bit resolution ADC
943 const Double_t vFSR = 1; // Full scale input voltage range
944 // const Double_t dNe = 125; // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
945 const Double_t dNe = 125/1.3; // F-ALTRO max V. FEE: factor 4
946 const Double_t vA = .136e-6; // CSP output range: 0.136uV/e-
947 const Double_t rise = 50e-9; // rise time (10-90%) of the FastOR signal before shaping
949 const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
951 Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
953 TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
954 signalF.SetParameter( 0, vV );
955 signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
956 signalF.SetParameter( 2, rise );
958 for (Int_t iTime=0; iTime<nSamples; iTime++)
960 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
961 // probably plan an access to OCDB
963 timeSamples[iTime] = int((TMath::Power(2, reso) / vFSR) * signalF.Eval(iTime * kTimeBinWidth) + 0.5);
967 //____________________________________________________________________________
968 Bool_t AliEMCALDigitizer::Init()
970 // Makes all memory allocations
972 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
974 if ( emcalLoader == 0 ) {
975 Fatal("Init", "Could not obtain the AliEMCALLoader");
980 fLastEvent = fFirstEvent ;
983 fInput = fDigInput->GetNinputs() ;
987 fInputFileNames = new TString[fInput] ;
988 fEventNames = new TString[fInput] ;
989 fInputFileNames[0] = GetTitle() ;
990 fEventNames[0] = fEventFolderName.Data() ;
992 for (index = 1 ; index < fInput ; index++) {
993 fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
994 TString tempo = fDigInput->GetInputFolderName(index) ;
995 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
998 //Calibration instance
999 fCalibData = emcalLoader->CalibData();
1003 //____________________________________________________________________________
1004 void AliEMCALDigitizer::InitParameters()
1006 // Parameter initialization for digitizer
1008 // Get the parameters from the OCDB via the loader
1009 AliRunLoader *rl = AliRunLoader::Instance();
1010 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1011 AliEMCALSimParam * simParam = 0x0;
1012 if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1015 simParam = AliEMCALSimParam::GetInstance();
1016 AliWarning("Simulation Parameters not available in OCDB?");
1019 fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400; // electrons per GeV
1020 fGainFluctuations = simParam->GetGainFluctuations() ; // 15.0;
1022 fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1023 if (fPinNoise < 0.0001 )
1024 Warning("InitParameters", "No noise added\n") ;
1025 fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
1026 fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1027 fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1028 fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1029 fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1031 // These defaults are normally not used.
1032 // Values are read from calibration database instead
1033 fADCchannelEC = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
1034 fADCpedestalEC = 0.0 ; // GeV
1035 fADCchannelECDecal = 1.0; // No decalibration by default
1036 fTimeChannel = 0.0; // No time calibration by default
1037 fTimeChannelDecal = 0.0; // No time decalibration by default
1039 fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
1041 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",
1042 fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1046 //__________________________________________________________________
1047 void AliEMCALDigitizer::Print1(Option_t * option)
1048 { // 19-nov-04 - just for convenience
1050 PrintDigits(option);
1053 //__________________________________________________________________
1054 void AliEMCALDigitizer::Print (Option_t * ) const
1056 // Print Digitizer's parameters
1057 printf("Print: \n------------------- %s -------------", GetName() ) ;
1058 if( strcmp(fEventFolderName.Data(), "") != 0 ){
1059 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1063 nStreams = GetNInputStreams() ;
1070 for (index = 0 ; index < nStreams ; index++) {
1071 TString tempo(fEventNames[index]) ;
1073 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1074 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1075 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1077 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1078 if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1079 fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
1080 printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1084 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1086 if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1087 else printf("\nNULL LOADER");
1089 printf("\nWith following parameters:\n") ;
1090 printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1091 printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
1092 printf("---------------------------------------------------\n") ;
1095 printf("Print: AliEMCALDigitizer not initialized") ;
1098 //__________________________________________________________________
1099 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1101 //utility method for printing digit information
1103 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1105 TClonesArray * digits = emcalLoader->Digits() ;
1106 TClonesArray * sdigits = emcalLoader->SDigits() ;
1108 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1109 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1111 if(strstr(option,"all")){
1113 AliEMCALDigit * digit;
1114 printf("\nEMC digits (with primaries):\n") ;
1115 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1117 for (index = 0 ; index < digits->GetEntries() ; index++) {
1118 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1120 printf("\n%6d %8f %6.5e %4d %2d : ",
1121 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1123 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1124 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1131 else printf("NULL LOADER, cannot print\n");
1134 //__________________________________________________________________
1135 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1137 // Calculates the time signal generated by noise
1138 //printf("Time noise %e\n",fTimeNoise);
1139 return gRandom->Rndm() * fTimeNoise;
1142 //__________________________________________________________________
1143 void AliEMCALDigitizer::Unload()
1145 // Unloads the SDigits and Digits
1149 for(i = 1 ; i < fInput ; i++){
1150 TString tempo(fEventNames[i]) ;
1152 if ((rl = AliRunLoader::GetRunLoader(tempo)))
1153 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1155 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1156 if(emcalLoader)emcalLoader->UnloadDigits() ;
1157 else AliFatal("Did not get the loader");
1160 //_________________________________________________________________________________________
1161 void AliEMCALDigitizer::WriteDigits()
1162 { // Makes TreeD in the output file.
1163 // Check if branch already exists:
1164 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1165 // already existing branches.
1166 // else creates branch with Digits, named "EMCAL", title "...",
1167 // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1168 // and names of files, from which digits are made.
1170 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1173 const TClonesArray * digits = emcalLoader->Digits() ;
1174 TTree * treeD = emcalLoader->TreeD();
1176 emcalLoader->MakeDigitsContainer();
1177 treeD = emcalLoader->TreeD();
1180 // -- create Digits branch
1181 Int_t bufferSize = 32000 ;
1182 TBranch * digitsBranch = 0;
1183 if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1184 digitsBranch->SetAddress(&digits);
1185 AliWarning("Digits Branch already exists. Not all digits will be visible");
1188 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1189 //digitsBranch->SetTitle(fEventFolderName);
1193 emcalLoader->WriteDigits("OVERWRITE");
1194 emcalLoader->WriteDigitizer("OVERWRITE");
1200 else AliFatal("Loader not available");
1203 //__________________________________________________________________
1204 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1205 { // overloaded method
1206 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1209 TTree* treeD = emcalLoader->TreeD();
1212 emcalLoader->MakeDigitsContainer();
1213 treeD = emcalLoader->TreeD();
1216 // -- create Digits branch
1217 Int_t bufferSize = 32000;
1219 if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1221 triggerBranch->SetAddress(&digits);
1225 treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1230 else AliFatal("Loader not available");