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"
70 #include "AliCaloCalibPedestal.h"
74 Double_t HeavisideTheta(Double_t x)
78 if (x > 0.) signal = 1.;
83 Double_t AnalogFastORFunction(Double_t *x, Double_t *par)
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));
109 ClassImp(AliEMCALDigitizer)
112 //____________________________________________________________________________
113 AliEMCALDigitizer::AliEMCALDigitizer()
114 : AliDigitizer("",""),
119 fInputFileNames(0x0),
122 fMeanPhotonElectron(0),
123 fGainFluctuations(0),
127 fTimeResolutionPar0(0),
128 fTimeResolutionPar1(0),
131 fADCchannelECDecal(0),
133 fTimeChannelDecal(0),
135 fEventFolderName(""),
143 fDigInput = 0 ; // We work in the standalone mode
146 //____________________________________________________________________________
147 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
148 : AliDigitizer("EMCALDigitizer", alirunFileName),
149 fDefaultInit(kFALSE),
156 fMeanPhotonElectron(0),
157 fGainFluctuations(0),
161 fTimeResolutionPar0(0),
162 fTimeResolutionPar1(0),
165 fADCchannelECDecal(0),
167 fTimeChannelDecal(0),
169 fEventFolderName(eventFolderName),
178 fDigInput = 0 ; // We work in the standalone mode
181 //____________________________________________________________________________
182 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
183 : AliDigitizer(d.GetName(),d.GetTitle()),
184 fDefaultInit(d.fDefaultInit),
185 fDigitsInRun(d.fDigitsInRun),
188 fInputFileNames(d.fInputFileNames),
189 fEventNames(d.fEventNames),
190 fDigitThreshold(d.fDigitThreshold),
191 fMeanPhotonElectron(d.fMeanPhotonElectron),
192 fGainFluctuations(d.fGainFluctuations),
193 fPinNoise(d.fPinNoise),
194 fTimeNoise(d.fTimeNoise),
195 fTimeDelay(d.fTimeDelay),
196 fTimeResolutionPar0(d.fTimeResolutionPar0),
197 fTimeResolutionPar1(d.fTimeResolutionPar1),
198 fADCchannelEC(d.fADCchannelEC),
199 fADCpedestalEC(d.fADCpedestalEC),
200 fADCchannelECDecal(d.fADCchannelECDecal),
201 fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
203 fEventFolderName(d.fEventFolderName),
204 fFirstEvent(d.fFirstEvent),
205 fLastEvent(d.fLastEvent),
206 fCalibData(d.fCalibData),
207 fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
212 //____________________________________________________________________________
213 AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
214 : AliDigitizer(rd,"EMCALDigitizer"),
215 fDefaultInit(kFALSE),
222 fMeanPhotonElectron(0),
223 fGainFluctuations(0),
227 fTimeResolutionPar0(0.),
228 fTimeResolutionPar1(0.),
231 fADCchannelECDecal(0),
233 fTimeChannelDecal(0),
241 // ctor Init() is called by RunDigitizer
243 fEventFolderName = fDigInput->GetInputFolderName(0) ;
244 SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
248 //____________________________________________________________________________
249 AliEMCALDigitizer::~AliEMCALDigitizer()
252 delete [] fInputFileNames ;
253 delete [] fEventNames ;
254 if (fSDigitizer) delete fSDigitizer;
257 //____________________________________________________________________________
258 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"));
277 Int_t readEvent = event ;
278 // fDigInput is data member from AliDigitizer
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())) ;
283 rl->GetEvent(readEvent);
285 TClonesArray * digits = emcalLoader->Digits() ;
286 digits->Delete() ; //correct way to clear array when memory is
287 //allocated by objects stored in array
290 AliEMCALGeometry *geom = 0;
291 if (!rl->GetAliRun())
293 AliFatal("Could not get AliRun from runLoader");
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()));
304 digits->Expand(nEMC) ;
306 // RS create a digitizer on the fly
307 if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
308 fSDigitizer->SetEventRange(0, -1) ;
310 //take all the inputs to add together and load the SDigits
311 TObjArray * sdigArray = new TObjArray(fInput) ;
312 sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
314 Int_t embed = kFALSE;
315 for(i = 1 ; i < fInput ; i++)
317 TString tempo(fEventNames[i]) ;
320 AliRunLoader * rl2 = AliRunLoader::GetRunLoader(tempo) ;
323 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
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()) ;
329 // rl2->LoadDigits();
330 rl2->GetEvent(readEvent);
331 if(rl2->GetDetectorLoader("EMCAL"))
333 AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
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)
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;
352 else AliFatal("EMCAL Loader of second event not available!");
355 else Fatal("Digitize", "Loader of second input not found");
358 //Find the first tower with signal
359 Int_t nextSig = nEMC + 1 ;
360 TClonesArray * sdigits ;
361 for(i = 0 ; i < fInput ; i++){
362 if(i > 0 && embed) continue;
363 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
366 if (sdigits->GetEntriesFast() )
368 AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
371 Int_t curNext = sd->GetId() ;
372 if(curNext < nextSig)
374 AliDebug(1,Form("input %i : #sdigits %i \n",
375 i, sdigits->GetEntriesFast()));
377 }//First entry is not NULL
380 AliDebug(1, "NULL sdigit pointer");
383 }//There is at least one entry
386 AliDebug(1, "NULL sdigits array");
389 }// SDigits array exists
392 AliDebug(1,"Null sdigit pointer");
396 AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
398 TArrayI index(fInput) ;
399 index.Reset() ; //Set all indexes to zero
401 AliEMCALDigit * digit ;
402 AliEMCALDigit * curSDigit ;
404 //Put Noise contribution, smear time and energy
405 Float_t timeResolution = 0;
406 for(absID = 0; absID < nEMC; absID++)
407 { // Nov 30, 2006 by PAI; was from 1 to nEMC
409 if (IsDead(absID)) continue; // Don't instantiate dead digits
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
424 // Calculate time as time of the largest digit
425 Float_t time = digit->GetTime() ;
426 Float_t aTime= digit->GetAmplitude() ;
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));
435 Int_t sDigitEntries = sdtclarr->GetEntriesFast();
437 if(sDigitEntries > index[i] )
438 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
442 //May be several digits will contribute from the same input
443 while(curSDigit && (curSDigit->GetId() == absID))
445 //Shift primary to separate primaries belonging different inputs
446 Int_t primaryoffset ;
448 primaryoffset = fDigInput->GetMask(i) ;
451 curSDigit->ShiftPrimary(primaryoffset) ;
453 if(curSDigit->GetAmplitude()>aTime)
455 aTime = curSDigit->GetAmplitude();
456 time = curSDigit->GetTime();
459 *digit = *digit + *curSDigit ; //adds amplitudes of each digit
463 if( sDigitEntries > index[i] )
464 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
469 }// loop over merging sources
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
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 ;
480 //calculate and set time
481 digit->SetTime(time) ;
483 //Find next signal module
485 for(i = 0 ; i < nInputs ; i++){
486 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
489 Int_t curNext = nextSig ;
490 if(sdigits->GetEntriesFast() > index[i])
492 AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
495 curNext = tmpdigit->GetId() ;
498 if(curNext < nextSig) nextSig = curNext ;
504 // add the noise now, no need for embedded events with real data
506 energy += gRandom->Gaus(0., fPinNoise) ;
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)) ;
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));
521 digit->SetTime(digit->GetTime()+fTimeDelay) ;
523 }// Digit pointer exists
524 else AliDebug(1,"Digit pointer is null");
526 } // for(absID = 0; absID < nEMC; absID++)
528 //Embed simulated amplitude (and time?) to real data digits
531 for(Int_t input = 1; input<fInput; input++)
533 TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
536 AliDebug(1,"No sdigits to merge\n");
540 for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
542 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
545 digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
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);
555 Float_t e0 = digit ->GetAmplitude();
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 ));
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);
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()));
578 //JLK is it better to call Clear() here?
579 delete sdigArray ; //We should not delete its contents
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
587 for(i = 0 ; i < nEMC ; i++)
589 digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
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) ;
605 Int_t ndigits = digits->GetEntriesFast() ;
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.
611 for (i = 0 ; i < ndigits ; i++)
613 digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
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());
628 else AliFatal("EMCALLoader is NULL!") ;
632 //_____________________________________________________________________
633 void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
636 // Returns digitized value of the energy in a cell absId
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
643 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
647 AliFatal("Did not get geometry from EMCALLoader");
657 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
659 Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
660 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
664 fADCpedestalEC = fCalibData->GetADCpedestal (iSupMod,ieta,iphi);
665 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
666 fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
667 fTimeChannel = fCalibData->GetTimeChannel (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
668 fTimeChannelDecal = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
671 //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
672 energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal ;
673 time += fTimeChannel-fTimeChannelDecal;
675 if(energy > fNADCEC )
679 //_____________________________________________________________________
680 void AliEMCALDigitizer::Decalibrate(AliEMCALDigit *digit)
682 // Decalibrate, used in Trigger digits
684 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
688 AliFatal("Did not get geometry from EMCALLoader");
692 Int_t absId = digit->GetId();
700 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
702 if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
703 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
707 fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
708 Float_t factor = 1./(fADCchannelEC/0.0162);
710 *digit = *digit * factor;
714 //_____________________________________________________________________
715 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
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
722 // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
723 // or with time out of window
726 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
730 AliFatal("Did not get geometry from EMCALLoader");
740 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
741 if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
742 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
746 fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
747 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
748 fTimeChannel = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
751 adc = adc * fADCchannelEC - fADCpedestalEC;
752 time -= fTimeChannel;
757 //____________________________________________________________________________
758 void AliEMCALDigitizer::Digitize(Option_t *option)
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)
766 if (!fInit) { // to prevent overwrite existing file
767 Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
771 if (strstr(option,"print")) {
777 if(strstr(option,"tim"))
778 gBenchmark->Start("EMCALDigitizer");
780 AliRunLoader *rl = AliRunLoader::Instance();
781 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
784 AliFatal("Did not get the Loader");
788 if (fLastEvent == -1) {
789 fLastEvent = rl->GetNumberOfEvents() - 1 ;
792 fLastEvent = fFirstEvent ; // what is this ??
794 nEvents = fLastEvent - fFirstEvent + 1;
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);
802 rl->LoadSDigits("EMCAL");
803 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
805 rl->GetEvent(ievent);
807 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
812 //-------------------------------------
815 Digits2FastOR(digitsTMP, digitsTRG);
817 WriteDigits(digitsTRG);
819 (emcalLoader->TreeD())->Fill();
821 emcalLoader->WriteDigits( "OVERWRITE");
825 digitsTRG ->Delete();
826 digitsTMP ->Delete();
828 //-------------------------------------
830 if(strstr(option,"deb"))
832 if(strstr(option,"table")) gObjectTable->Print();
834 //increment the total number of Digits per run
835 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
840 if(strstr(option,"tim")){
841 gBenchmark->Stop("EMCALDigitizer");
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)) ;
849 //__________________________________________________________________
850 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
852 // Assign a smeared time to the digit, from observed distribution in LED system (?).
857 res = TMath::Sqrt(fTimeResolutionPar0 +
858 fTimeResolutionPar1/(energy*energy) );
860 // parametrization above is for ns. Convert to seconds:
864 //____________________________________________________________________________
865 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
867 // FEE digits afterburner to produce TRG digits
868 // we are only interested in the FEE digit deposited energy
869 // to be converted later into a voltage value
871 // push the FEE digit to its associated FastOR (numbered from 0:95)
872 // TRU is in charge of summing module digits
874 AliRunLoader *runLoader = AliRunLoader::Instance();
876 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
877 if (!emcalLoader) AliFatal("Did not get the Loader");
879 const AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
881 // build FOR from simulated digits
882 // and xfer to the corresponding TRU input (mapping)
884 TClonesArray* sdigits = emcalLoader->SDigits();
886 TClonesArray *digits = (TClonesArray*)sdigits->Clone();
888 AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
890 TIter NextDigit(digits);
891 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
893 if (IsDead(digit)) continue;
897 Int_t id = digit->GetId();
900 if (geom && geom->GetFastORIndexFromCellIndex(id , trgid))
902 AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
904 AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
908 new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
909 d = (AliEMCALDigit*)digitsTMP->At(trgid);
919 if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
921 Int_t nSamples = geom->GetNTotalTRU();
922 Int_t *timeSamples = new Int_t[nSamples];
924 NextDigit = TIter(digitsTMP);
925 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
929 Int_t id = digit->GetId();
930 Float_t time = 50.e-9;
932 Double_t depositedEnergy = 0.;
933 for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
935 if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
937 // FIXME: Check digit time!
938 if (depositedEnergy) {
939 depositedEnergy += gRandom->Gaus(0., .08);
940 DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
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));
947 new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
949 if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
954 delete [] timeSamples;
956 if (digits && digits->GetEntriesFast()) digits->Delete();
959 //____________________________________________________________________________
960 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
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
971 const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
973 Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
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 );
980 for (Int_t iTime=0; iTime<nSamples; iTime++)
982 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
983 // probably plan an access to OCDB
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;
989 AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
990 timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
995 //____________________________________________________________________________
996 Bool_t AliEMCALDigitizer::Init()
998 // Makes all memory allocations
1000 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1002 if ( emcalLoader == 0 ) {
1003 Fatal("Init", "Could not obtain the AliEMCALLoader");
1008 fLastEvent = fFirstEvent ;
1011 fInput = fDigInput->GetNinputs() ;
1015 fInputFileNames = new TString[fInput] ;
1016 fEventNames = new TString[fInput] ;
1017 fInputFileNames[0] = GetTitle() ;
1018 fEventNames[0] = fEventFolderName.Data() ;
1020 for (index = 1 ; index < fInput ; index++) {
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
1026 //Calibration instance
1027 fCalibData = emcalLoader->CalibData();
1031 //____________________________________________________________________________
1032 void AliEMCALDigitizer::InitParameters()
1034 // Parameter initialization for digitizer
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();
1043 simParam = AliEMCALSimParam::GetInstance();
1044 AliWarning("Simulation Parameters not available in OCDB?");
1047 fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400; // electrons per GeV
1048 fGainFluctuations = simParam->GetGainFluctuations() ; // 15.0;
1050 fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1051 if (fPinNoise < 0.0001 )
1052 Warning("InitParameters", "No noise added\n") ;
1053 fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
1054 fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1055 fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1056 fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1057 fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1059 // These defaults are normally not used.
1060 // Values are read from calibration database instead
1061 fADCchannelEC = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
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
1067 fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
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));
1074 //__________________________________________________________________
1075 void AliEMCALDigitizer::Print1(Option_t * option)
1076 { // 19-nov-04 - just for convenience
1078 PrintDigits(option);
1081 //__________________________________________________________________
1082 void AliEMCALDigitizer::Print (Option_t * ) const
1084 // Print Digitizer's parameters
1085 printf("Print: \n------------------- %s -------------", GetName() ) ;
1086 if( strcmp(fEventFolderName.Data(), "") != 0 ){
1087 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1091 nStreams = GetNInputStreams() ;
1098 for (index = 0 ; index < nStreams ; index++) {
1099 TString tempo(fEventNames[index]) ;
1101 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1102 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1103 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
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()) ;
1112 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1114 if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1115 else printf("\nNULL LOADER");
1117 printf("\nWith following parameters:\n") ;
1118 printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1119 printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
1120 printf("---------------------------------------------------\n") ;
1123 printf("Print: AliEMCALDigitizer not initialized") ;
1126 //__________________________________________________________________
1127 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1129 //utility method for printing digit information
1131 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1133 TClonesArray * digits = emcalLoader->Digits() ;
1134 TClonesArray * sdigits = emcalLoader->SDigits() ;
1136 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1137 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1139 if(strstr(option,"all")){
1141 AliEMCALDigit * digit;
1142 printf("\nEMC digits (with primaries):\n") ;
1143 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1145 for (index = 0 ; index < digits->GetEntries() ; index++) {
1146 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1148 printf("\n%6d %8f %6.5e %4d %2d : ",
1149 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1151 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1152 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1159 else printf("NULL LOADER, cannot print\n");
1162 //__________________________________________________________________
1163 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1165 // Calculates the time signal generated by noise
1166 //printf("Time noise %e\n",fTimeNoise);
1167 return gRandom->Rndm() * fTimeNoise;
1170 //__________________________________________________________________
1171 void AliEMCALDigitizer::Unload()
1173 // Unloads the SDigits and Digits
1177 for(i = 1 ; i < fInput ; i++){
1178 TString tempo(fEventNames[i]) ;
1180 if ((rl = AliRunLoader::GetRunLoader(tempo)))
1181 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1183 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1184 if(emcalLoader)emcalLoader->UnloadDigits() ;
1185 else AliFatal("Did not get the loader");
1188 //_________________________________________________________________________________________
1189 void AliEMCALDigitizer::WriteDigits()
1190 { // Makes TreeD in the output file.
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.
1198 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1201 const TClonesArray * digits = emcalLoader->Digits() ;
1202 TTree * treeD = emcalLoader->TreeD();
1204 emcalLoader->MakeDigitsContainer();
1205 treeD = emcalLoader->TreeD();
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");
1216 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1217 //digitsBranch->SetTitle(fEventFolderName);
1221 emcalLoader->WriteDigits("OVERWRITE");
1222 emcalLoader->WriteDigitizer("OVERWRITE");
1228 else AliFatal("Loader not available");
1231 //__________________________________________________________________
1232 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1233 { // overloaded method
1234 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1237 TTree* treeD = emcalLoader->TreeD();
1240 emcalLoader->MakeDigitsContainer();
1241 treeD = emcalLoader->TreeD();
1244 // -- create Digits branch
1245 Int_t bufferSize = 32000;
1247 if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1249 triggerBranch->SetAddress(&digits);
1253 treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1258 else AliFatal("Loader not available");
1261 //__________________________________________________________________
1262 Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1264 AliRunLoader *runLoader = AliRunLoader::Instance();
1265 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1266 if (!emcalLoader) AliFatal("Did not get the Loader");
1268 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1270 AliWarning("Could not access pedestal data! No dead channel removal applied");
1275 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1276 if (!geom) AliFatal("Did not get geometry from EMCALLoader");
1278 Int_t absId = digit->GetId();
1286 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1288 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1289 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1291 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1293 if (channelStatus == AliCaloCalibPedestal::kDead)
1300 //__________________________________________________________________
1301 Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1303 AliRunLoader *runLoader = AliRunLoader::Instance();
1304 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1305 if (!emcalLoader) AliFatal("Did not get the Loader");
1307 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1309 AliWarning("Could not access pedestal data! No dead channel removal applied");
1314 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1315 if (!geom) AliFatal("Did not get geometry from EMCALLoader");
1324 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1326 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1327 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1329 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1331 if (channelStatus == AliCaloCalibPedestal::kDead)