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),
125 fTimeResolutionPar0(0),
126 fTimeResolutionPar1(0),
129 fADCchannelECDecal(0),
131 fTimeChannelDecal(0),
133 fEventFolderName(""),
141 fDigInput = 0 ; // We work in the standalone mode
144 //____________________________________________________________________________
145 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
146 : AliDigitizer("EMCALDigitizer", alirunFileName),
147 fDefaultInit(kFALSE),
154 fMeanPhotonElectron(0),
158 fTimeResolutionPar0(0),
159 fTimeResolutionPar1(0),
162 fADCchannelECDecal(0),
164 fTimeChannelDecal(0),
166 fEventFolderName(eventFolderName),
175 fDigInput = 0 ; // We work in the standalone mode
178 //____________________________________________________________________________
179 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
180 : AliDigitizer(d.GetName(),d.GetTitle()),
181 fDefaultInit(d.fDefaultInit),
182 fDigitsInRun(d.fDigitsInRun),
185 fInputFileNames(d.fInputFileNames),
186 fEventNames(d.fEventNames),
187 fDigitThreshold(d.fDigitThreshold),
188 fMeanPhotonElectron(d.fMeanPhotonElectron),
189 fPinNoise(d.fPinNoise),
190 fTimeNoise(d.fTimeNoise),
191 fTimeDelay(d.fTimeDelay),
192 fTimeResolutionPar0(d.fTimeResolutionPar0),
193 fTimeResolutionPar1(d.fTimeResolutionPar1),
194 fADCchannelEC(d.fADCchannelEC),
195 fADCpedestalEC(d.fADCpedestalEC),
196 fADCchannelECDecal(d.fADCchannelECDecal),
197 fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
199 fEventFolderName(d.fEventFolderName),
200 fFirstEvent(d.fFirstEvent),
201 fLastEvent(d.fLastEvent),
202 fCalibData(d.fCalibData),
203 fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
208 //____________________________________________________________________________
209 AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
210 : AliDigitizer(rd,"EMCALDigitizer"),
211 fDefaultInit(kFALSE),
218 fMeanPhotonElectron(0),
222 fTimeResolutionPar0(0.),
223 fTimeResolutionPar1(0.),
226 fADCchannelECDecal(0),
228 fTimeChannelDecal(0),
236 // ctor Init() is called by RunDigitizer
238 fEventFolderName = fDigInput->GetInputFolderName(0) ;
239 SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
243 //____________________________________________________________________________
244 AliEMCALDigitizer::~AliEMCALDigitizer()
247 delete [] fInputFileNames ;
248 delete [] fEventNames ;
249 if (fSDigitizer) delete fSDigitizer;
252 //____________________________________________________________________________
253 void AliEMCALDigitizer::Digitize(Int_t event)
256 // Makes the digitization of the collected summable digits
257 // for this it first creates the array of all EMCAL modules
258 // filled with noise and after that adds contributions from
259 // SDigits. This design helps to avoid scanning over the
260 // list of digits to add contribution of any new SDigit.
263 // Note that SDigit energy info is stored as an amplitude, so we
264 // must call the Calibrate() method of the SDigitizer to convert it
265 // back to an energy in GeV before adding it to the Digit
267 static int nEMC=0; //max number of digits possible
268 AliRunLoader *rl = AliRunLoader::Instance();
269 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
272 Int_t readEvent = event ;
273 // fDigInput is data member from AliDigitizer
275 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
276 AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
277 readEvent, GetTitle(), fEventFolderName.Data())) ;
278 rl->GetEvent(readEvent);
280 TClonesArray * digits = emcalLoader->Digits() ;
281 digits->Delete() ; //correct way to clear array when memory is
282 //allocated by objects stored in array
285 AliEMCALGeometry *geom = 0;
286 if (!rl->GetAliRun()) {
287 AliFatal("Could not get AliRun from runLoader");
290 AliEMCAL * emcal = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
291 geom = emcal->GetGeometry();
292 nEMC = geom->GetNCells();
293 AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
298 digits->Expand(nEMC) ;
300 // RS create a digitizer on the fly
301 if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
302 fSDigitizer->SetEventRange(0, -1) ;
304 //take all the inputs to add together and load the SDigits
305 TObjArray * sdigArray = new TObjArray(fInput) ;
306 sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
308 Int_t embed = kFALSE;
309 for(i = 1 ; i < fInput ; i++){
310 TString tempo(fEventNames[i]) ;
313 AliRunLoader * rl2 = AliRunLoader::GetRunLoader(tempo) ;
316 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
319 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
320 Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
322 // rl2->LoadDigits();
323 rl2->GetEvent(readEvent);
324 if(rl2->GetDetectorLoader("EMCAL"))
326 AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
329 //Merge 2 simulated sdigits
330 if(emcalLoader2->SDigits()){
331 TClonesArray* sdigits2 = emcalLoader2->SDigits();
332 sdigArray->AddAt(sdigits2, i) ;
333 // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
334 // do not smear energy of embedded, do not add noise to any sdigits
335 if(sdigits2->GetEntriesFast()>0){
336 //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
337 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
338 if(digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded){
345 else AliFatal("EMCAL Loader of second event not available!");
348 else Fatal("Digitize", "Loader of second input not found");
351 //Find the first tower with signal
352 Int_t nextSig = nEMC + 1 ;
353 TClonesArray * sdigits ;
354 for(i = 0 ; i < fInput ; i++){
355 if(i > 0 && embed) continue;
356 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
358 if (sdigits->GetEntriesFast() ){
359 AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
361 Int_t curNext = sd->GetId() ;
362 if(curNext < nextSig)
364 AliDebug(1,Form("input %i : #sdigits %i \n",
365 i, sdigits->GetEntriesFast()));
367 }//First entry is not NULL
369 Info("Digitize", "NULL sdigit pointer");
372 }//There is at least one entry
374 Info("Digitize", "NULL sdigits array");
377 }// SDigits array exists
379 Info("Digitizer","Null sdigit pointer");
383 AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
385 TArrayI index(fInput) ;
386 index.Reset() ; //Set all indexes to zero
388 AliEMCALDigit * digit ;
389 AliEMCALDigit * curSDigit ;
391 //Put Noise contribution, smear time and energy
392 Float_t timeResolution = 0;
393 for(absID = 0; absID < nEMC; absID++){ // Nov 30, 2006 by PAI; was from 1 to nEMC
396 // amplitude set to zero, noise will be added later
397 Float_t noiseTime = 0.;
398 if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
399 new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
400 //look if we have to add signal?
401 digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
407 // Calculate time as time of the largest digit
408 Float_t time = digit->GetTime() ;
409 Float_t aTime= digit->GetAmplitude() ;
412 Int_t nInputs = fInput;
413 if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
414 for(i = 0; i< nInputs ; i++){ //loop over (possible) merge sources
415 TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
417 Int_t sDigitEntries = sdtclarr->GetEntriesFast();
419 if(sDigitEntries > index[i] )
420 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
424 //May be several digits will contribute from the same input
425 while(curSDigit && (curSDigit->GetId() == absID)){
426 //Shift primary to separate primaries belonging different inputs
427 Int_t primaryoffset ;
429 primaryoffset = fDigInput->GetMask(i) ;
432 curSDigit->ShiftPrimary(primaryoffset) ;
434 if(curSDigit->GetAmplitude()>aTime) {
435 aTime = curSDigit->GetAmplitude();
436 time = curSDigit->GetTime();
439 *digit = *digit + *curSDigit ; //adds amplitudes of each digit
443 if( sDigitEntries > index[i] )
444 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
449 }// loop over merging sources
452 //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
453 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
455 // add fluctuations for photo-electron creation
456 energy *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
458 //calculate and set time
459 digit->SetTime(time) ;
461 //Find next signal module
463 for(i = 0 ; i < nInputs ; i++){
464 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
467 Int_t curNext = nextSig ;
468 if(sdigits->GetEntriesFast() > index[i])
470 AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
473 curNext = tmpdigit->GetId() ;
476 if(curNext < nextSig) nextSig = curNext ;
482 // add the noise now, no need for embedded events with real data
484 energy += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
488 //Now digitize the energy according to the fSDigitizer method,
489 //which merely converts the energy to an integer. Later we will
490 //check that the stored value matches our allowed dynamic ranges
491 digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
493 //Set time resolution with final energy
494 timeResolution = GetTimeResolution(energy);
495 digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
496 AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
497 absID, energy, nextSig));
499 digit->SetTime(digit->GetTime()+fTimeDelay) ;
501 }// Digit pointer exists
503 Info("Digitizer","Digit pointer is null");
505 } // for(absID = 0; absID < nEMC; absID++)
510 //Embed simulated amplitude (and time?) to real data digits
512 for(Int_t input = 1; input<fInput; input++){
513 TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
515 AliDebug(1,"No sdigits to merge\n");
518 for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++){
519 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
521 digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
524 // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
525 // and multiply to get a big int.
526 Float_t time2 = digit2->GetTime();
527 Float_t e2 = digit2->GetAmplitude();
528 CalibrateADCTime(e2,time2,digit2->GetId());
529 Float_t amp2 = fSDigitizer->Digitize(e2);
531 Float_t e0 = digit ->GetAmplitude();
533 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",
534 digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
535 digit2->GetId(),amp2, digit2->GetType(), time2 ));
537 // Sum amplitudes, change digit amplitude
538 digit->SetAmplitude( digit->GetAmplitude() + amp2);
539 // Rough assumption, assign time of the largest of the 2 digits,
540 // data or signal to the final digit.
541 if(amp2 > digit->GetAmplitude()) digit->SetTime(time2);
542 //Mark the new digit as embedded
543 digit->SetType(AliEMCALDigit::kEmbedded);
545 if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
546 AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
547 digit->GetId(), digit->GetAmplitude(), digit->GetType()));
554 //JLK is it better to call Clear() here?
555 delete sdigArray ; //We should not delete its contents
557 //remove digits below thresholds
558 // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
559 // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
560 // before merge in the same loop real data digits if available
563 for(i = 0 ; i < nEMC ; i++){
564 digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
567 //Then get the energy in GeV units.
568 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
569 //Then digitize using the calibration constants of the ocdb
570 Float_t ampADC = energy;
571 DigitizeEnergyTime(ampADC, time, digit->GetId()) ;
572 if(ampADC < fDigitThreshold)
573 digits->RemoveAt(i) ;
580 Int_t ndigits = digits->GetEntriesFast() ;
583 //After we have done the summing and digitizing to create the
584 //digits, now we want to calibrate the resulting amplitude to match
585 //the dynamic range of our real data.
586 for (i = 0 ; i < ndigits ; i++) {
587 digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
589 digit->SetIndexInList(i) ;
590 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
591 time = digit->GetTime();
592 Float_t ampADC = energy;
593 DigitizeEnergyTime(ampADC, time, digit->GetId());
594 digit->SetAmplitude(ampADC) ;
595 digit->SetTime(time);
596 // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
601 else AliFatal("EMCALLoader is NULL!") ;
605 //_____________________________________________________________________
606 void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
609 // Returns digitized value of the energy in a cell absId
610 // using the calibration constants stored in the OCDB
611 // or default values if no CalibData object is found.
612 // This effectively converts everything to match the dynamic range
613 // of the real data we will collect
616 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
619 AliFatal("Did not get geometry from EMCALLoader");
629 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
631 Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
632 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
635 fADCpedestalEC = fCalibData->GetADCpedestal (iSupMod,ieta,iphi);
636 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
637 fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
638 fTimeChannel = fCalibData->GetTimeChannel (iSupMod,ieta,iphi);
639 fTimeChannelDecal = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
642 //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
643 energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal ;
644 time += fTimeChannel-fTimeChannelDecal;
646 if(energy > fNADCEC )
650 //_____________________________________________________________________
651 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
653 // Returns the energy in a cell absId with a given adc value
654 // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
655 // Used in case of embedding, transform ADC counts from real event into energy
656 // so that we can add the energy of the simulated sdigits which are in energy
658 // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
659 // or with time out of window
662 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
665 AliFatal("Did not get geometry from EMCALLoader");
675 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
677 Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
678 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
681 fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
682 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
683 fTimeChannel = fCalibData->GetTimeChannel(iSupMod,ieta,iphi);
686 adc = adc * fADCchannelEC - fADCpedestalEC;
687 time -= fTimeChannel;
693 //____________________________________________________________________________
694 void AliEMCALDigitizer::Digitize(Option_t *option)
696 // Steering method to process digitization for events
697 // in the range from fFirstEvent to fLastEvent.
698 // This range is optionally set by SetEventRange().
699 // if fLastEvent=-1, then process events until the end.
700 // by default fLastEvent = fFirstEvent (process only one event)
702 if (!fInit) { // to prevent overwrite existing file
703 Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
707 if (strstr(option,"print")) {
713 if(strstr(option,"tim"))
714 gBenchmark->Start("EMCALDigitizer");
716 AliRunLoader *rl = AliRunLoader::Instance();
717 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
720 AliFatal("Did not get the Loader");
724 if (fLastEvent == -1) {
725 fLastEvent = rl->GetNumberOfEvents() - 1 ;
728 fLastEvent = fFirstEvent ; // what is this ??
730 nEvents = fLastEvent - fFirstEvent + 1;
733 TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit", 32*96);
734 TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", 32*96);
736 rl->LoadSDigits("EMCAL");
737 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
739 rl->GetEvent(ievent);
741 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
746 //-------------------------------------
749 Digits2FastOR(digitsTMP, digitsTRG);
751 WriteDigits(digitsTRG);
753 (emcalLoader->TreeD())->Fill();
755 emcalLoader->WriteDigits( "OVERWRITE");
759 digitsTRG ->Delete();
760 digitsTMP ->Delete();
762 //-------------------------------------
764 if(strstr(option,"deb"))
766 if(strstr(option,"table")) gObjectTable->Print();
768 //increment the total number of Digits per run
769 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
774 if(strstr(option,"tim")){
775 gBenchmark->Stop("EMCALDigitizer");
776 AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event",
777 gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
781 //__________________________________________________________________
782 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
787 res = TMath::Sqrt(fTimeResolutionPar0 +
788 fTimeResolutionPar1/(energy*energy) );
790 // parametrization above is for ns. Convert to seconds:
794 //____________________________________________________________________________
795 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
797 // FEE digits afterburner to produce TRG digits
798 // we are only interested in the FEE digit deposited energy
799 // to be converted later into a voltage value
801 // push the FEE digit to its associated FastOR (numbered from 0:95)
802 // TRU is in charge of summing module digits
804 AliRunLoader *runLoader = AliRunLoader::Instance();
806 AliRun* run = runLoader->GetAliRun();
808 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
810 AliFatal("Did not get the Loader");
813 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(run->GetDetector("EMCAL"));
815 AliEMCALGeometry* geom = emcal->GetGeometry();
817 // build FOR from simulated digits
818 // and xfer to the corresponding TRU input (mapping)
820 TClonesArray* digits = emcalLoader->Digits();
822 TIter NextDigit(digits);
823 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
825 Int_t id = digit->GetId();
828 if (geom && geom->GetFastORIndexFromCellIndex(id , trgid))
830 AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
832 AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
836 new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
837 d = (AliEMCALDigit*)digitsTMP->At(trgid);
847 if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
850 Int_t *timeSamples = new Int_t[nSamples];
852 NextDigit = TIter(digitsTMP);
853 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
857 Int_t id = digit->GetId();
858 Float_t time = digit->GetTime();
860 Double_t depositedEnergy = 0.;
861 for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
863 if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
865 // FIXME: Check digit time!
868 DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
870 for (Int_t j=0;j<nSamples;j++)
872 timeSamples[j] = ((j << 12) & 0xFF000) | (timeSamples[j] & 0xFFF);
875 new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
880 delete [] timeSamples;
882 else AliFatal("Could not get AliEMCAL");
887 //____________________________________________________________________________
888 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
892 const Int_t reso = 11; // 11-bit resolution ADC
893 const Double_t vFSR = 1; // Full scale input voltage range
894 const Double_t dNe = 125; // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
895 const Double_t vA = .136e-6; // CSP output range: 0.136uV/e-
896 const Double_t rise = 40e-9; // rise time (10-90%) of the FastOR signal before shaping
898 const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
900 Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
902 TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
903 signalF.SetParameter( 0, vV );
904 signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
905 signalF.SetParameter( 2, rise );
907 for (Int_t iTime=0; iTime<nSamples; iTime++)
909 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
910 // probably plan an access to OCDB
912 timeSamples[iTime] = int((TMath::Power(2, reso) / vFSR) * signalF.Eval(iTime * kTimeBinWidth) + 0.5);
916 //____________________________________________________________________________
917 Bool_t AliEMCALDigitizer::Init()
919 // Makes all memory allocations
921 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
923 if ( emcalLoader == 0 ) {
924 Fatal("Init", "Could not obtain the AliEMCALLoader");
929 fLastEvent = fFirstEvent ;
932 fInput = fDigInput->GetNinputs() ;
936 fInputFileNames = new TString[fInput] ;
937 fEventNames = new TString[fInput] ;
938 fInputFileNames[0] = GetTitle() ;
939 fEventNames[0] = fEventFolderName.Data() ;
941 for (index = 1 ; index < fInput ; index++) {
942 fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
943 TString tempo = fDigInput->GetInputFolderName(index) ;
944 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
947 //Calibration instance
948 fCalibData = emcalLoader->CalibData();
952 //____________________________________________________________________________
953 void AliEMCALDigitizer::InitParameters()
955 // Parameter initialization for digitizer
957 // Get the parameters from the OCDB via the loader
958 AliRunLoader *rl = AliRunLoader::Instance();
959 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
960 AliEMCALSimParam * simParam = 0x0;
961 if(emcalLoader) simParam = emcalLoader->SimulationParameters();
964 simParam = AliEMCALSimParam::GetInstance();
965 AliWarning("Simulation Parameters not available in OCDB?");
968 fMeanPhotonElectron = simParam->GetMeanPhotonElectron();//4400; // electrons per GeV
969 fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
970 if (fPinNoise < 0.0001 )
971 Warning("InitParameters", "No noise added\n") ;
972 fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
973 fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
974 fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
975 fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
976 fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
978 // These defaults are normally not used.
979 // Values are read from calibration database instead
980 fADCchannelEC = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
981 fADCpedestalEC = 0.0 ; // GeV
982 fADCchannelECDecal = 1.0; // No decalibration by default
983 fTimeChannel = 0.0; // No time calibration by default
984 fTimeChannelDecal = 0.0; // No time decalibration by default
986 fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
988 AliDebug(2,Form("Mean Photon Electron %d, Noise: APD %f, Time %f; Digit Threshold %d,Time Resolution Par0 %g Par1 %g,NADCEC %d",
989 fMeanPhotonElectron,fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
993 //__________________________________________________________________
994 void AliEMCALDigitizer::Print1(Option_t * option)
995 { // 19-nov-04 - just for convenience
1000 //__________________________________________________________________
1001 void AliEMCALDigitizer::Print (Option_t * ) const
1003 // Print Digitizer's parameters
1004 printf("Print: \n------------------- %s -------------", GetName() ) ;
1005 if( strcmp(fEventFolderName.Data(), "") != 0 ){
1006 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1010 nStreams = GetNInputStreams() ;
1017 for (index = 0 ; index < nStreams ; index++) {
1018 TString tempo(fEventNames[index]) ;
1020 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1021 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1022 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1024 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1025 if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1026 fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
1027 printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1031 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1033 if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1034 else printf("\nNULL LOADER");
1036 printf("\nWith following parameters:\n") ;
1037 printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1038 printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
1039 printf("---------------------------------------------------\n") ;
1042 printf("Print: AliEMCALDigitizer not initialized") ;
1045 //__________________________________________________________________
1046 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1048 //utility method for printing digit information
1050 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1052 TClonesArray * digits = emcalLoader->Digits() ;
1053 TClonesArray * sdigits = emcalLoader->SDigits() ;
1055 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1056 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1058 if(strstr(option,"all")){
1060 AliEMCALDigit * digit;
1061 printf("\nEMC digits (with primaries):\n") ;
1062 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1064 for (index = 0 ; index < digits->GetEntries() ; index++) {
1065 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1067 printf("\n%6d %8f %6.5e %4d %2d : ",
1068 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1070 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1071 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1078 else printf("NULL LOADER, cannot print\n");
1081 //__________________________________________________________________
1082 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1084 // Calculates the time signal generated by noise
1085 //printf("Time noise %e\n",fTimeNoise);
1086 return gRandom->Rndm() * fTimeNoise;
1089 //__________________________________________________________________
1090 void AliEMCALDigitizer::Unload()
1092 // Unloads the SDigits and Digits
1096 for(i = 1 ; i < fInput ; i++){
1097 TString tempo(fEventNames[i]) ;
1099 if ((rl = AliRunLoader::GetRunLoader(tempo)))
1100 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1102 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1103 if(emcalLoader)emcalLoader->UnloadDigits() ;
1104 else AliFatal("Did not get the loader");
1107 //_________________________________________________________________________________________
1108 void AliEMCALDigitizer::WriteDigits()
1109 { // Makes TreeD in the output file.
1110 // Check if branch already exists:
1111 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1112 // already existing branches.
1113 // else creates branch with Digits, named "EMCAL", title "...",
1114 // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1115 // and names of files, from which digits are made.
1117 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1120 const TClonesArray * digits = emcalLoader->Digits() ;
1121 TTree * treeD = emcalLoader->TreeD();
1123 emcalLoader->MakeDigitsContainer();
1124 treeD = emcalLoader->TreeD();
1127 // -- create Digits branch
1128 Int_t bufferSize = 32000 ;
1129 TBranch * digitsBranch = 0;
1130 if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1131 digitsBranch->SetAddress(&digits);
1132 AliWarning("Digits Branch already exists. Not all digits will be visible");
1135 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1136 //digitsBranch->SetTitle(fEventFolderName);
1140 emcalLoader->WriteDigits("OVERWRITE");
1141 emcalLoader->WriteDigitizer("OVERWRITE");
1147 else AliFatal("Loader not available");
1150 //__________________________________________________________________
1151 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1152 { // overloaded method
1153 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1156 TTree* treeD = emcalLoader->TreeD();
1159 emcalLoader->MakeDigitsContainer();
1160 treeD = emcalLoader->TreeD();
1163 // -- create Digits branch
1164 Int_t bufferSize = 32000;
1166 if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1168 triggerBranch->SetAddress(&digits);
1172 treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1177 else AliFatal("Loader not available");