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
272 AliRunLoader *rl = AliRunLoader::Instance();
273 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
277 AliFatal("EMCALLoader is NULL!") ;
278 return; // not needed but in case coverity complains ...
281 Int_t readEvent = event ;
282 if (fDigInput) // fDigInput is data member from AliDigitizer
283 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
284 AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
285 readEvent, GetTitle(), fEventFolderName.Data())) ;
286 rl->GetEvent(readEvent);
288 TClonesArray * digits = emcalLoader->Digits() ;
289 digits->Delete() ; //correct way to clear array when memory is
290 //allocated by objects stored in array
293 if (!rl->GetAliRun())
295 AliFatal("Could not get AliRun from runLoader");
296 return; // not needed but in case coverity complains ...
299 AliEMCAL * emcal = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
300 AliEMCALGeometry *geom = emcal->GetGeometry();
301 static int nEMC = geom->GetNCells();//max number of digits possible
302 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 //-------------------------------------------------------
311 //take all the inputs to add together and load the SDigits
312 TObjArray * sdigArray = new TObjArray(fInput) ;
313 sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
316 Int_t embed = kFALSE;
317 for(i = 1 ; i < fInput ; i++)
319 TString tempo(fEventNames[i]) ;
322 AliRunLoader * rl2 = AliRunLoader::GetRunLoader(tempo) ;
325 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
329 AliFatal("Run Loader of second event not available!");
330 return; // not needed but in case coverity complains ...
334 readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
336 Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ;
339 // rl2->LoadDigits();
340 rl2->GetEvent(readEvent);
342 if(!rl2->GetDetectorLoader("EMCAL"))
344 AliFatal("Loader of second input not found");
345 return; // not needed but in case coverity complains ...
348 AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
352 AliFatal("EMCAL Loader of second event not available!");
353 return; // not needed but in case coverity complains ...
356 //Merge 2 simulated sdigits
357 if(!emcalLoader2->SDigits()) continue;
359 TClonesArray* sdigits2 = emcalLoader2->SDigits();
360 sdigArray->AddAt(sdigits2, i) ;
362 // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
363 // do not smear energy of embedded, do not add noise to any sdigits
364 if( sdigits2->GetEntriesFast() <= 0 ) continue;
366 //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
367 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
368 if( digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded ) embed = kTRUE;
372 //--------------------------------
373 //Find the first tower with signal
374 Int_t nextSig = nEMC + 1 ;
375 TClonesArray * sdigits ;
376 for(i = 0 ; i < fInput ; i++)
378 if(i > 0 && embed) continue;
380 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
383 AliDebug(1,"Null sdigit pointer");
387 if (sdigits->GetEntriesFast() <= 0 )
389 AliDebug(1, "No sdigits entries");
393 AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
396 AliDebug(1, "NULL sdigit pointer");
400 Int_t curNext = sd->GetId() ;
401 if(curNext < nextSig)
403 AliDebug(1,Form("input %i : #sdigits %i \n",i, sdigits->GetEntriesFast()));
407 AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
409 TArrayI index(fInput) ;
410 index.Reset() ; //Set all indexes to zero
412 AliEMCALDigit * digit ;
413 AliEMCALDigit * curSDigit ;
415 //---------------------------------------------
416 //Put Noise contribution, smear time and energy
417 Float_t timeResolution = 0;
419 for(absID = 0; absID < nEMC; absID++)
423 // amplitude set to zero, noise will be added later
424 Float_t noiseTime = 0.;
425 if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
426 new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
427 //look if we have to add signal?
428 digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
432 AliDebug(1,"Digit pointer is null");
438 // Calculate time as time of the largest digit
439 Float_t time = digit->GetTime() ;
440 Float_t aTime= digit->GetAmplitude() ;
443 Int_t nInputs = fInput;
444 if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
445 for(i = 0; i< nInputs ; i++)
447 //loop over (possible) merge sources
448 TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
451 Int_t sDigitEntries = sdtclarr->GetEntriesFast();
453 if(sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
456 //May be several digits will contribute from the same input
457 while(curSDigit && (curSDigit->GetId() == absID))
459 //Shift primary to separate primaries belonging different inputs
460 Int_t primaryoffset = i ;
461 if(fDigInput) primaryoffset = fDigInput->GetMask(i) ;
462 curSDigit->ShiftPrimary(primaryoffset) ;
464 if(curSDigit->GetAmplitude()>aTime)
466 aTime = curSDigit->GetAmplitude();
467 time = curSDigit->GetTime();
470 *digit = *digit + *curSDigit ; //adds amplitudes of each digit
474 if( sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
478 }// loop over merging sources
480 //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
481 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
483 // add fluctuations for photo-electron creation
484 // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
485 Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
486 energy *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
488 //calculate and set time
489 digit->SetTime(time) ;
491 //Find next signal module
493 for(i = 0 ; i < nInputs ; i++)
495 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
499 Int_t curNext = nextSig ;
500 if(sdigits->GetEntriesFast() > index[i])
502 AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
503 if ( tmpdigit ) curNext = tmpdigit->GetId() ;
506 if(curNext < nextSig) nextSig = curNext ;
512 // add the noise now, no need for embedded events with real data
514 energy += gRandom->Gaus(0., fPinNoise) ;
517 //Now digitize the energy according to the fSDigitizer method,
518 //which merely converts the energy to an integer. Later we will
519 //check that the stored value matches our allowed dynamic ranges
520 digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
522 //Set time resolution with final energy
523 timeResolution = GetTimeResolution(energy);
524 digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
525 AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
526 absID, energy, nextSig));
528 digit->SetTime(digit->GetTime()+fTimeDelay) ;
530 } // for(absID = 0; absID < nEMC; absID++)
532 //---------------------------------------------------------
533 //Embed simulated amplitude (and time?) to real data digits
536 for(Int_t input = 1; input<fInput; input++)
538 TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
541 AliDebug(1,"No sdigits to merge\n");
545 for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
547 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
548 if ( !digit2 ) continue;
550 digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
551 if ( !digit ) continue;
553 // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
554 // and multiply to get a big int.
555 Float_t time2 = digit2->GetTime();
556 Float_t e2 = digit2->GetAmplitude();
557 CalibrateADCTime(e2,time2,digit2->GetId());
558 Float_t amp2 = fSDigitizer->Digitize(e2);
560 Float_t e0 = digit ->GetAmplitude();
562 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",
563 digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
564 digit2->GetId(),amp2, digit2->GetType(), time2 ));
566 // Sum amplitudes, change digit amplitude
567 digit->SetAmplitude( digit->GetAmplitude() + amp2);
569 // Rough assumption, assign time of the largest of the 2 digits,
570 // data or signal to the final digit.
571 if(amp2 > digit->GetAmplitude()) digit->SetTime(time2);
573 //Mark the new digit as embedded
574 digit->SetType(AliEMCALDigit::kEmbedded);
576 if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
577 AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
578 digit->GetId(), digit->GetAmplitude(), digit->GetType()));
583 //JLK is it better to call Clear() here?
584 delete sdigArray ; //We should not delete its contents
586 //------------------------------
587 //remove digits below thresholds
588 // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
589 // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
590 // before merge in the same loop real data digits if available
593 for(i = 0 ; i < nEMC ; i++)
595 digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
596 if ( !digit ) continue;
598 //Then get the energy in GeV units.
599 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
601 //Then digitize using the calibration constants of the ocdb
602 Float_t ampADC = energy;
603 DigitizeEnergyTime(ampADC, time, digit->GetId()) ;
605 if(ampADC < fDigitThreshold || IsDead(digit->GetId()))
606 digits->RemoveAt(i) ;
612 Int_t ndigits = digits->GetEntriesFast() ;
614 //---------------------------------------------------------------
616 //After we have done the summing and digitizing to create the
617 //digits, now we want to calibrate the resulting amplitude to match
618 //the dynamic range of our real data.
619 for (i = 0 ; i < ndigits ; i++)
621 digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
622 if( !digit ) continue ;
624 digit->SetIndexInList(i) ;
626 time = digit->GetTime();
627 digit->SetTime(time);
629 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
631 Float_t ampADC = energy;
632 DigitizeEnergyTime(ampADC, time, digit->GetId());
634 digit->SetAmplitude(ampADC) ;
635 // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
640 //_____________________________________________________________________
641 void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
644 // Returns digitized value of the energy in a cell absId
645 // using the calibration constants stored in the OCDB
646 // or default values if no CalibData object is found.
647 // This effectively converts everything to match the dynamic range
648 // of the real data we will collect
651 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
655 AliFatal("Did not get geometry from EMCALLoader");
665 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
667 Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
668 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
672 fADCpedestalEC = fCalibData->GetADCpedestal (iSupMod,ieta,iphi);
673 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
674 fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
675 fTimeChannel = fCalibData->GetTimeChannel (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
676 fTimeChannelDecal = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
679 //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
680 energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal ;
681 time += fTimeChannel-fTimeChannelDecal;
683 if ( energy > fNADCEC ) energy = fNADCEC ;
686 //_____________________________________________________________________
687 void AliEMCALDigitizer::DecalibrateTrigger(AliEMCALDigit *digit)
689 // Decalibrate, used in Trigger digits
691 if ( !fCalibData ) return ;
694 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
698 AliFatal("Did not get geometry from EMCALLoader");
702 // Recover the digit location in row-column-SM
703 Int_t absId = digit->GetId();
711 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
712 if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
714 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
717 Float_t adcChannel = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
718 Float_t adcChannelOnline = fCalibData->GetADCchannelOnline(iSupMod,ieta,iphi);
719 //printf("calib param for (SM%d,col%d,row%d): %1.4f - online param: %1.4f \n",iSupMod,ieta,iphi,adcChannel,adcChannelOnline);
720 Float_t factor = 1./(adcChannel/adcChannelOnline);
722 *digit = *digit * factor;
726 //_____________________________________________________________________
727 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
729 // Returns the energy in a cell absId with a given adc value
730 // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
731 // Used in case of embedding, transform ADC counts from real event into energy
732 // so that we can add the energy of the simulated sdigits which are in energy
734 // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
735 // or with time out of window
738 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
742 AliFatal("Did not get geometry from EMCALLoader");
752 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
753 if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
754 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
758 fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
759 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
760 fTimeChannel = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
763 adc = adc * fADCchannelEC - fADCpedestalEC;
764 time -= fTimeChannel;
769 //____________________________________________________________________________
770 void AliEMCALDigitizer::Digitize(Option_t *option)
772 // Steering method to process digitization for events
773 // in the range from fFirstEvent to fLastEvent.
774 // This range is optionally set by SetEventRange().
775 // if fLastEvent=-1, then process events until the end.
776 // by default fLastEvent = fFirstEvent (process only one event)
779 { // to prevent overwrite existing file
780 Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
784 if (strstr(option,"print"))
790 if(strstr(option,"tim"))
791 gBenchmark->Start("EMCALDigitizer");
793 AliRunLoader *rl = AliRunLoader::Instance();
795 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
798 AliFatal("Did not get the Loader");
802 if (fLastEvent == -1)
803 fLastEvent = rl->GetNumberOfEvents() - 1 ;
805 fLastEvent = fFirstEvent ; // what is this ??
807 Int_t nEvents = fLastEvent - fFirstEvent + 1;
810 AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
813 AliFatal("Did not get the AliEMCAL pointer");
817 AliEMCALGeometry *geom = emcal->GetGeometry();
820 AliFatal("Geometry pointer null");
821 return; // fix for coverity
824 const Int_t nTRU = geom->GetNTotalTRU();
825 TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit", nTRU*96);
826 TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
828 rl->LoadSDigits("EMCAL");
829 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++)
831 rl->GetEvent(ievent);
833 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
838 //-------------------------------------
840 Digits2FastOR(digitsTMP, digitsTRG);
842 WriteDigits(digitsTRG);
844 (emcalLoader->TreeD())->Fill();
846 emcalLoader->WriteDigits( "OVERWRITE");
850 digitsTRG ->Delete();
851 digitsTMP ->Delete();
853 //-------------------------------------
855 if(strstr(option,"deb"))
857 if(strstr(option,"table")) gObjectTable->Print();
859 //increment the total number of Digits per run
860 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
863 if(strstr(option,"tim"))
865 gBenchmark->Stop("EMCALDigitizer");
866 Float_t cputime = gBenchmark->GetCpuTime("EMCALDigitizer");
867 Float_t avcputime = cputime;
868 if(nEvents==0) avcputime = 0 ;
869 AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
873 //__________________________________________________________________
874 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
876 // Assign a smeared time to the digit, from observed distribution in LED system (?).
881 res = TMath::Sqrt(fTimeResolutionPar0 +
882 fTimeResolutionPar1/(energy*energy) );
884 // parametrization above is for ns. Convert to seconds:
888 //____________________________________________________________________________
889 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
891 // FEE digits afterburner to produce TRG digits
892 // we are only interested in the FEE digit deposited energy
893 // to be converted later into a voltage value
895 // push the FEE digit to its associated FastOR (numbered from 0:95)
896 // TRU is in charge of summing module digits
898 AliRunLoader *runLoader = AliRunLoader::Instance();
900 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
901 if (!emcalLoader) AliFatal("Did not get the Loader");
903 const AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
906 AliFatal("Geometry pointer null");
907 return; // fix for coverity
910 // build FOR from simulated digits
911 // and xfer to the corresponding TRU input (mapping)
913 TClonesArray* sdigits = emcalLoader->SDigits();
915 TClonesArray *digits = (TClonesArray*)sdigits->Clone();
917 AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
919 TIter NextDigit(digits);
920 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
922 if (IsDead(digit)) continue;
924 DecalibrateTrigger(digit);
926 Int_t id = digit->GetId();
929 if (geom->GetFastORIndexFromCellIndex(id , trgid))
931 AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
933 AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
937 new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
938 d = (AliEMCALDigit*)digitsTMP->At(trgid);
948 if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
950 Int_t nSamples = geom->GetNTotalTRU();
951 Int_t *timeSamples = new Int_t[nSamples];
953 NextDigit = TIter(digitsTMP);
954 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
958 Int_t id = digit->GetId();
959 Float_t time = 50.e-9;
961 Double_t depositedEnergy = 0.;
962 for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
964 if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
966 // FIXME: Check digit time!
967 if (depositedEnergy) {
968 depositedEnergy += gRandom->Gaus(0., .08);
969 DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
971 for (Int_t j=0;j<nSamples;j++) {
972 if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
973 timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
976 new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
978 if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
983 delete [] timeSamples;
985 if (digits && digits->GetEntriesFast()) digits->Delete();
988 //____________________________________________________________________________
989 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
993 const Int_t reso = 12; // 11-bit resolution ADC
994 const Double_t vFSR = 2.; // Full scale input voltage range 2V (p-p)
995 //const Double_t dNe = 125; // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
996 const Double_t dNe = 125/1.3; // F-ALTRO max V. FEE: factor 4
997 const Double_t vA = .136e-6; // CSP output range: 0.136uV/e-
998 const Double_t rise = 50e-9; // rise time (10-90%) of the FastOR signal before shaping
1000 const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
1002 Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
1004 TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
1005 signalF.SetParameter( 0, vV );
1006 signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
1007 signalF.SetParameter( 2, rise );
1009 for (Int_t iTime=0; iTime<nSamples; iTime++)
1011 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
1012 // probably plan an access to OCDB
1013 Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
1014 if (TMath::Abs(sig) > vFSR/2.) {
1015 AliError("Signal overflow!");
1016 timeSamples[iTime] = (1 << reso) - 1;
1018 AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
1019 timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
1024 //____________________________________________________________________________
1025 Bool_t AliEMCALDigitizer::Init()
1027 // Makes all memory allocations
1029 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1031 if ( emcalLoader == 0 ) {
1032 Fatal("Init", "Could not obtain the AliEMCALLoader");
1037 fLastEvent = fFirstEvent ;
1040 fInput = fDigInput->GetNinputs() ;
1044 fInputFileNames = new TString[fInput] ;
1045 fEventNames = new TString[fInput] ;
1046 fInputFileNames[0] = GetTitle() ;
1047 fEventNames[0] = fEventFolderName.Data() ;
1049 for (index = 1 ; index < fInput ; index++) {
1050 fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1051 TString tempo = fDigInput->GetInputFolderName(index) ;
1052 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
1055 //Calibration instance
1056 fCalibData = emcalLoader->CalibData();
1060 //____________________________________________________________________________
1061 void AliEMCALDigitizer::InitParameters()
1063 // Parameter initialization for digitizer
1065 // Get the parameters from the OCDB via the loader
1066 AliRunLoader *rl = AliRunLoader::Instance();
1067 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1068 AliEMCALSimParam * simParam = 0x0;
1069 if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1072 simParam = AliEMCALSimParam::GetInstance();
1073 AliWarning("Simulation Parameters not available in OCDB?");
1076 fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400; // electrons per GeV
1077 fGainFluctuations = simParam->GetGainFluctuations() ; // 15.0;
1079 fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1080 if (fPinNoise < 0.0001 )
1081 Warning("InitParameters", "No noise added\n") ;
1082 fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
1083 fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1084 fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1085 fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1086 fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1088 // These defaults are normally not used.
1089 // Values are read from calibration database instead
1090 fADCchannelEC = 0.0162; // Nominal value set online for most of the channels, MIP peak sitting at ~266./16/1024
1091 fADCpedestalEC = 0.0 ; // GeV
1092 fADCchannelECDecal = 1.0; // No decalibration by default
1093 fTimeChannel = 0.0; // No time calibration by default
1094 fTimeChannelDecal = 0.0; // No time decalibration by default
1096 fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
1098 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",
1099 fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1103 //__________________________________________________________________
1104 void AliEMCALDigitizer::Print1(Option_t * option)
1105 { // 19-nov-04 - just for convenience
1107 PrintDigits(option);
1110 //__________________________________________________________________
1111 void AliEMCALDigitizer::Print (Option_t * ) const
1113 // Print Digitizer's parameters
1114 printf("Print: \n------------------- %s -------------", GetName() ) ;
1115 if( strcmp(fEventFolderName.Data(), "") != 0 ){
1116 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1120 nStreams = GetNInputStreams() ;
1127 for (index = 0 ; index < nStreams ; index++) {
1128 TString tempo(fEventNames[index]) ;
1130 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1131 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1132 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1134 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1135 if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1136 fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
1137 printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1141 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1143 if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1144 else printf("\nNULL LOADER");
1146 printf("\nWith following parameters:\n") ;
1147 printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1148 printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
1149 printf("---------------------------------------------------\n") ;
1152 printf("Print: AliEMCALDigitizer not initialized") ;
1155 //__________________________________________________________________
1156 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1158 //utility method for printing digit information
1160 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1162 TClonesArray * digits = emcalLoader->Digits() ;
1163 TClonesArray * sdigits = emcalLoader->SDigits() ;
1165 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1166 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1168 if(strstr(option,"all")){
1170 AliEMCALDigit * digit;
1171 printf("\nEMC digits (with primaries):\n") ;
1172 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1174 for (index = 0 ; index < digits->GetEntries() ; index++) {
1175 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1177 printf("\n%6d %8f %6.5e %4d %2d : ",
1178 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1180 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1181 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1188 else printf("NULL LOADER, cannot print\n");
1191 //__________________________________________________________________
1192 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1194 // Calculates the time signal generated by noise
1195 //printf("Time noise %e\n",fTimeNoise);
1196 return gRandom->Rndm() * fTimeNoise;
1199 //__________________________________________________________________
1200 void AliEMCALDigitizer::Unload()
1202 // Unloads the SDigits and Digits
1206 for(i = 1 ; i < fInput ; i++){
1207 TString tempo(fEventNames[i]) ;
1209 if ((rl = AliRunLoader::GetRunLoader(tempo)))
1210 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1212 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1213 if(emcalLoader)emcalLoader->UnloadDigits() ;
1214 else AliFatal("Did not get the loader");
1217 //_________________________________________________________________________________________
1218 void AliEMCALDigitizer::WriteDigits()
1219 { // Makes TreeD in the output file.
1220 // Check if branch already exists:
1221 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1222 // already existing branches.
1223 // else creates branch with Digits, named "EMCAL", title "...",
1224 // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1225 // and names of files, from which digits are made.
1227 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1230 const TClonesArray * digits = emcalLoader->Digits() ;
1231 TTree * treeD = emcalLoader->TreeD();
1233 emcalLoader->MakeDigitsContainer();
1234 treeD = emcalLoader->TreeD();
1237 // -- create Digits branch
1238 Int_t bufferSize = 32000 ;
1239 TBranch * digitsBranch = 0;
1240 if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1241 digitsBranch->SetAddress(&digits);
1242 AliWarning("Digits Branch already exists. Not all digits will be visible");
1245 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1246 //digitsBranch->SetTitle(fEventFolderName);
1250 emcalLoader->WriteDigits("OVERWRITE");
1251 emcalLoader->WriteDigitizer("OVERWRITE");
1257 else AliFatal("Loader not available");
1260 //__________________________________________________________________
1261 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1263 // overloaded method
1264 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1267 TTree* treeD = emcalLoader->TreeD();
1270 emcalLoader->MakeDigitsContainer();
1271 treeD = emcalLoader->TreeD();
1274 // -- create Digits branch
1275 Int_t bufferSize = 32000;
1277 if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1279 triggerBranch->SetAddress(&digits);
1283 treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1288 else AliFatal("Loader not available");
1291 //__________________________________________________________________
1292 Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1294 // Check if cell is defined as dead, so that it is not included
1296 Int_t absId = digit->GetId();
1298 return IsDead(absId);
1303 //__________________________________________________________________
1304 Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1306 // Check if cell absID is defined as dead, so that it is not included
1308 AliRunLoader *runLoader = AliRunLoader::Instance();
1309 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1312 AliFatal("Did not get the Loader");
1316 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1319 AliWarning("Could not access pedestal data! No dead channel removal applied");
1324 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1327 AliFatal("Did not get geometry from EMCALLoader");
1328 return kTRUE; //coverity
1338 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1340 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1341 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1343 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1345 if (channelStatus == AliCaloCalibPedestal::kDead)