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++)
420 { // Nov 30, 2006 by PAI; was from 1 to nEMC
422 if (IsDead(absID)) continue; // Don't instantiate dead digits
426 // amplitude set to zero, noise will be added later
427 Float_t noiseTime = 0.;
428 if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
429 new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
430 //look if we have to add signal?
431 digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
435 AliDebug(1,"Digit pointer is null");
441 // Calculate time as time of the largest digit
442 Float_t time = digit->GetTime() ;
443 Float_t aTime= digit->GetAmplitude() ;
446 Int_t nInputs = fInput;
447 if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
448 for(i = 0; i< nInputs ; i++)
450 //loop over (possible) merge sources
451 TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
454 Int_t sDigitEntries = sdtclarr->GetEntriesFast();
456 if(sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
459 //May be several digits will contribute from the same input
460 while(curSDigit && (curSDigit->GetId() == absID))
462 //Shift primary to separate primaries belonging different inputs
463 Int_t primaryoffset = i ;
464 if(fDigInput) primaryoffset = fDigInput->GetMask(i) ;
465 curSDigit->ShiftPrimary(primaryoffset) ;
467 if(curSDigit->GetAmplitude()>aTime)
469 aTime = curSDigit->GetAmplitude();
470 time = curSDigit->GetTime();
473 *digit = *digit + *curSDigit ; //adds amplitudes of each digit
477 if( sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
481 }// loop over merging sources
483 //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
484 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
486 // add fluctuations for photo-electron creation
487 // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
488 Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
489 energy *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
491 //calculate and set time
492 digit->SetTime(time) ;
494 //Find next signal module
496 for(i = 0 ; i < nInputs ; i++)
498 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
502 Int_t curNext = nextSig ;
503 if(sdigits->GetEntriesFast() > index[i])
505 AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
506 if ( tmpdigit ) curNext = tmpdigit->GetId() ;
509 if(curNext < nextSig) nextSig = curNext ;
515 // add the noise now, no need for embedded events with real data
517 energy += gRandom->Gaus(0., fPinNoise) ;
520 //Now digitize the energy according to the fSDigitizer method,
521 //which merely converts the energy to an integer. Later we will
522 //check that the stored value matches our allowed dynamic ranges
523 digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
525 //Set time resolution with final energy
526 timeResolution = GetTimeResolution(energy);
527 digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
528 AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
529 absID, energy, nextSig));
531 digit->SetTime(digit->GetTime()+fTimeDelay) ;
533 } // for(absID = 0; absID < nEMC; absID++)
535 //---------------------------------------------------------
536 //Embed simulated amplitude (and time?) to real data digits
539 for(Int_t input = 1; input<fInput; input++)
541 TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
544 AliDebug(1,"No sdigits to merge\n");
548 for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
550 AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
551 if ( !digit2 ) continue;
553 digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
554 if ( !digit ) continue;
556 // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
557 // and multiply to get a big int.
558 Float_t time2 = digit2->GetTime();
559 Float_t e2 = digit2->GetAmplitude();
560 CalibrateADCTime(e2,time2,digit2->GetId());
561 Float_t amp2 = fSDigitizer->Digitize(e2);
563 Float_t e0 = digit ->GetAmplitude();
565 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",
566 digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
567 digit2->GetId(),amp2, digit2->GetType(), time2 ));
569 // Sum amplitudes, change digit amplitude
570 digit->SetAmplitude( digit->GetAmplitude() + amp2);
572 // Rough assumption, assign time of the largest of the 2 digits,
573 // data or signal to the final digit.
574 if(amp2 > digit->GetAmplitude()) digit->SetTime(time2);
576 //Mark the new digit as embedded
577 digit->SetType(AliEMCALDigit::kEmbedded);
579 if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
580 AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
581 digit->GetId(), digit->GetAmplitude(), digit->GetType()));
586 //JLK is it better to call Clear() here?
587 delete sdigArray ; //We should not delete its contents
589 //------------------------------
590 //remove digits below thresholds
591 // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
592 // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
593 // before merge in the same loop real data digits if available
596 for(i = 0 ; i < nEMC ; i++)
598 digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
599 if ( !digit ) continue;
601 //Then get the energy in GeV units.
602 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
604 //Then digitize using the calibration constants of the ocdb
605 Float_t ampADC = energy;
606 DigitizeEnergyTime(ampADC, time, digit->GetId()) ;
608 if(ampADC < fDigitThreshold)
609 digits->RemoveAt(i) ;
615 Int_t ndigits = digits->GetEntriesFast() ;
617 //---------------------------------------------------------------
619 //After we have done the summing and digitizing to create the
620 //digits, now we want to calibrate the resulting amplitude to match
621 //the dynamic range of our real data.
622 for (i = 0 ; i < ndigits ; i++)
624 digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ;
625 if( !digit ) continue ;
627 digit->SetIndexInList(i) ;
629 time = digit->GetTime();
630 digit->SetTime(time);
632 energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
634 Float_t ampADC = energy;
635 DigitizeEnergyTime(ampADC, time, digit->GetId());
637 digit->SetAmplitude(ampADC) ;
638 // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
643 //_____________________________________________________________________
644 void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
647 // Returns digitized value of the energy in a cell absId
648 // using the calibration constants stored in the OCDB
649 // or default values if no CalibData object is found.
650 // This effectively converts everything to match the dynamic range
651 // of the real data we will collect
654 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
658 AliFatal("Did not get geometry from EMCALLoader");
668 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
670 Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
671 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
675 fADCpedestalEC = fCalibData->GetADCpedestal (iSupMod,ieta,iphi);
676 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
677 fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
678 fTimeChannel = fCalibData->GetTimeChannel (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
679 fTimeChannelDecal = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);
682 //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
683 energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal ;
684 time += fTimeChannel-fTimeChannelDecal;
686 if ( energy > fNADCEC ) energy = fNADCEC ;
689 //_____________________________________________________________________
690 void AliEMCALDigitizer::DecalibrateTrigger(AliEMCALDigit *digit)
692 // Decalibrate, used in Trigger digits
694 if ( !fCalibData ) return ;
697 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
701 AliFatal("Did not get geometry from EMCALLoader");
705 // Recover the digit location in row-column-SM
706 Int_t absId = digit->GetId();
714 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
715 if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
717 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
720 Float_t adcChannel = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
721 Float_t adcChannelOnline = fCalibData->GetADCchannelOnline(iSupMod,ieta,iphi);
722 //printf("calib param for (SM%d,col%d,row%d): %1.4f - online param: %1.4f \n",iSupMod,ieta,iphi,adcChannel,adcChannelOnline);
723 Float_t factor = 1./(adcChannel/adcChannelOnline);
725 *digit = *digit * factor;
729 //_____________________________________________________________________
730 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
732 // Returns the energy in a cell absId with a given adc value
733 // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
734 // Used in case of embedding, transform ADC counts from real event into energy
735 // so that we can add the energy of the simulated sdigits which are in energy
737 // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
738 // or with time out of window
741 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
745 AliFatal("Did not get geometry from EMCALLoader");
755 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
756 if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
757 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
761 fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
762 fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
763 fTimeChannel = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
766 adc = adc * fADCchannelEC - fADCpedestalEC;
767 time -= fTimeChannel;
772 //____________________________________________________________________________
773 void AliEMCALDigitizer::Digitize(Option_t *option)
775 // Steering method to process digitization for events
776 // in the range from fFirstEvent to fLastEvent.
777 // This range is optionally set by SetEventRange().
778 // if fLastEvent=-1, then process events until the end.
779 // by default fLastEvent = fFirstEvent (process only one event)
782 { // to prevent overwrite existing file
783 Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
787 if (strstr(option,"print"))
793 if(strstr(option,"tim"))
794 gBenchmark->Start("EMCALDigitizer");
796 AliRunLoader *rl = AliRunLoader::Instance();
798 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
801 AliFatal("Did not get the Loader");
805 if (fLastEvent == -1)
806 fLastEvent = rl->GetNumberOfEvents() - 1 ;
808 fLastEvent = fFirstEvent ; // what is this ??
810 Int_t nEvents = fLastEvent - fFirstEvent + 1;
813 AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
816 AliFatal("Did not get the AliEMCAL pointer");
820 AliEMCALGeometry *geom = emcal->GetGeometry();
823 AliFatal("Geometry pointer null");
824 return; // fix for coverity
827 const Int_t nTRU = geom->GetNTotalTRU();
828 TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit", nTRU*96);
829 TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
831 rl->LoadSDigits("EMCAL");
832 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++)
834 rl->GetEvent(ievent);
836 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
841 //-------------------------------------
843 Digits2FastOR(digitsTMP, digitsTRG);
845 WriteDigits(digitsTRG);
847 (emcalLoader->TreeD())->Fill();
849 emcalLoader->WriteDigits( "OVERWRITE");
853 digitsTRG ->Delete();
854 digitsTMP ->Delete();
856 //-------------------------------------
858 if(strstr(option,"deb"))
860 if(strstr(option,"table")) gObjectTable->Print();
862 //increment the total number of Digits per run
863 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
866 if(strstr(option,"tim"))
868 gBenchmark->Stop("EMCALDigitizer");
869 Float_t cputime = gBenchmark->GetCpuTime("EMCALDigitizer");
870 Float_t avcputime = cputime;
871 if(nEvents==0) avcputime = 0 ;
872 AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
876 //__________________________________________________________________
877 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
879 // Assign a smeared time to the digit, from observed distribution in LED system (?).
884 res = TMath::Sqrt(fTimeResolutionPar0 +
885 fTimeResolutionPar1/(energy*energy) );
887 // parametrization above is for ns. Convert to seconds:
891 //____________________________________________________________________________
892 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
894 // FEE digits afterburner to produce TRG digits
895 // we are only interested in the FEE digit deposited energy
896 // to be converted later into a voltage value
898 // push the FEE digit to its associated FastOR (numbered from 0:95)
899 // TRU is in charge of summing module digits
901 AliRunLoader *runLoader = AliRunLoader::Instance();
903 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
904 if (!emcalLoader) AliFatal("Did not get the Loader");
906 const AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
909 AliFatal("Geometry pointer null");
910 return; // fix for coverity
913 // build FOR from simulated digits
914 // and xfer to the corresponding TRU input (mapping)
916 TClonesArray* sdigits = emcalLoader->SDigits();
918 TClonesArray *digits = (TClonesArray*)sdigits->Clone();
920 AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
922 TIter NextDigit(digits);
923 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
925 if (IsDead(digit)) continue;
927 DecalibrateTrigger(digit);
929 Int_t id = digit->GetId();
932 if (geom->GetFastORIndexFromCellIndex(id , trgid))
934 AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
936 AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
940 new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
941 d = (AliEMCALDigit*)digitsTMP->At(trgid);
951 if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
953 Int_t nSamples = geom->GetNTotalTRU();
954 Int_t *timeSamples = new Int_t[nSamples];
956 NextDigit = TIter(digitsTMP);
957 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
961 Int_t id = digit->GetId();
962 Float_t time = 50.e-9;
964 Double_t depositedEnergy = 0.;
965 for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
967 if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
969 // FIXME: Check digit time!
970 if (depositedEnergy) {
971 depositedEnergy += gRandom->Gaus(0., .08);
972 DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
974 for (Int_t j=0;j<nSamples;j++) {
975 if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
976 timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
979 new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
981 if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
986 delete [] timeSamples;
988 if (digits && digits->GetEntriesFast()) digits->Delete();
991 //____________________________________________________________________________
992 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
996 const Int_t reso = 12; // 11-bit resolution ADC
997 const Double_t vFSR = 2.; // Full scale input voltage range 2V (p-p)
998 //const Double_t dNe = 125; // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
999 const Double_t dNe = 125/1.3; // F-ALTRO max V. FEE: factor 4
1000 const Double_t vA = .136e-6; // CSP output range: 0.136uV/e-
1001 const Double_t rise = 50e-9; // rise time (10-90%) of the FastOR signal before shaping
1003 const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
1005 Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
1007 TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
1008 signalF.SetParameter( 0, vV );
1009 signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
1010 signalF.SetParameter( 2, rise );
1012 for (Int_t iTime=0; iTime<nSamples; iTime++)
1014 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
1015 // probably plan an access to OCDB
1016 Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
1017 if (TMath::Abs(sig) > vFSR/2.) {
1018 AliError("Signal overflow!");
1019 timeSamples[iTime] = (1 << reso) - 1;
1021 AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
1022 timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
1027 //____________________________________________________________________________
1028 Bool_t AliEMCALDigitizer::Init()
1030 // Makes all memory allocations
1032 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1034 if ( emcalLoader == 0 ) {
1035 Fatal("Init", "Could not obtain the AliEMCALLoader");
1040 fLastEvent = fFirstEvent ;
1043 fInput = fDigInput->GetNinputs() ;
1047 fInputFileNames = new TString[fInput] ;
1048 fEventNames = new TString[fInput] ;
1049 fInputFileNames[0] = GetTitle() ;
1050 fEventNames[0] = fEventFolderName.Data() ;
1052 for (index = 1 ; index < fInput ; index++) {
1053 fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1054 TString tempo = fDigInput->GetInputFolderName(index) ;
1055 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
1058 //Calibration instance
1059 fCalibData = emcalLoader->CalibData();
1063 //____________________________________________________________________________
1064 void AliEMCALDigitizer::InitParameters()
1066 // Parameter initialization for digitizer
1068 // Get the parameters from the OCDB via the loader
1069 AliRunLoader *rl = AliRunLoader::Instance();
1070 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1071 AliEMCALSimParam * simParam = 0x0;
1072 if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1075 simParam = AliEMCALSimParam::GetInstance();
1076 AliWarning("Simulation Parameters not available in OCDB?");
1079 fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400; // electrons per GeV
1080 fGainFluctuations = simParam->GetGainFluctuations() ; // 15.0;
1082 fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1083 if (fPinNoise < 0.0001 )
1084 Warning("InitParameters", "No noise added\n") ;
1085 fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
1086 fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1087 fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1088 fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1089 fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1091 // These defaults are normally not used.
1092 // Values are read from calibration database instead
1093 fADCchannelEC = 0.0162; // Nominal value set online for most of the channels, MIP peak sitting at ~266./16/1024
1094 fADCpedestalEC = 0.0 ; // GeV
1095 fADCchannelECDecal = 1.0; // No decalibration by default
1096 fTimeChannel = 0.0; // No time calibration by default
1097 fTimeChannelDecal = 0.0; // No time decalibration by default
1099 fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
1101 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",
1102 fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1106 //__________________________________________________________________
1107 void AliEMCALDigitizer::Print1(Option_t * option)
1108 { // 19-nov-04 - just for convenience
1110 PrintDigits(option);
1113 //__________________________________________________________________
1114 void AliEMCALDigitizer::Print (Option_t * ) const
1116 // Print Digitizer's parameters
1117 printf("Print: \n------------------- %s -------------", GetName() ) ;
1118 if( strcmp(fEventFolderName.Data(), "") != 0 ){
1119 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1123 nStreams = GetNInputStreams() ;
1130 for (index = 0 ; index < nStreams ; index++) {
1131 TString tempo(fEventNames[index]) ;
1133 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1134 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1135 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1137 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1138 if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1139 fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
1140 printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1144 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1146 if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1147 else printf("\nNULL LOADER");
1149 printf("\nWith following parameters:\n") ;
1150 printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1151 printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
1152 printf("---------------------------------------------------\n") ;
1155 printf("Print: AliEMCALDigitizer not initialized") ;
1158 //__________________________________________________________________
1159 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1161 //utility method for printing digit information
1163 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1165 TClonesArray * digits = emcalLoader->Digits() ;
1166 TClonesArray * sdigits = emcalLoader->SDigits() ;
1168 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1169 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1171 if(strstr(option,"all")){
1173 AliEMCALDigit * digit;
1174 printf("\nEMC digits (with primaries):\n") ;
1175 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1177 for (index = 0 ; index < digits->GetEntries() ; index++) {
1178 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1180 printf("\n%6d %8f %6.5e %4d %2d : ",
1181 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1183 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1184 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1191 else printf("NULL LOADER, cannot print\n");
1194 //__________________________________________________________________
1195 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1197 // Calculates the time signal generated by noise
1198 //printf("Time noise %e\n",fTimeNoise);
1199 return gRandom->Rndm() * fTimeNoise;
1202 //__________________________________________________________________
1203 void AliEMCALDigitizer::Unload()
1205 // Unloads the SDigits and Digits
1209 for(i = 1 ; i < fInput ; i++){
1210 TString tempo(fEventNames[i]) ;
1212 if ((rl = AliRunLoader::GetRunLoader(tempo)))
1213 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1215 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1216 if(emcalLoader)emcalLoader->UnloadDigits() ;
1217 else AliFatal("Did not get the loader");
1220 //_________________________________________________________________________________________
1221 void AliEMCALDigitizer::WriteDigits()
1222 { // Makes TreeD in the output file.
1223 // Check if branch already exists:
1224 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1225 // already existing branches.
1226 // else creates branch with Digits, named "EMCAL", title "...",
1227 // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1228 // and names of files, from which digits are made.
1230 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1233 const TClonesArray * digits = emcalLoader->Digits() ;
1234 TTree * treeD = emcalLoader->TreeD();
1236 emcalLoader->MakeDigitsContainer();
1237 treeD = emcalLoader->TreeD();
1240 // -- create Digits branch
1241 Int_t bufferSize = 32000 ;
1242 TBranch * digitsBranch = 0;
1243 if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1244 digitsBranch->SetAddress(&digits);
1245 AliWarning("Digits Branch already exists. Not all digits will be visible");
1248 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1249 //digitsBranch->SetTitle(fEventFolderName);
1253 emcalLoader->WriteDigits("OVERWRITE");
1254 emcalLoader->WriteDigitizer("OVERWRITE");
1260 else AliFatal("Loader not available");
1263 //__________________________________________________________________
1264 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1266 // overloaded method
1267 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1270 TTree* treeD = emcalLoader->TreeD();
1273 emcalLoader->MakeDigitsContainer();
1274 treeD = emcalLoader->TreeD();
1277 // -- create Digits branch
1278 Int_t bufferSize = 32000;
1280 if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1282 triggerBranch->SetAddress(&digits);
1286 treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1291 else AliFatal("Loader not available");
1294 //__________________________________________________________________
1295 Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1297 // Check if cell is defined as dead, so that it is not included
1299 Int_t absId = digit->GetId();
1301 return IsDead(absId);
1306 //__________________________________________________________________
1307 Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1309 // Check if cell absID is defined as dead, so that it is not included
1311 AliRunLoader *runLoader = AliRunLoader::Instance();
1312 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1315 AliFatal("Did not get the Loader");
1319 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1322 AliWarning("Could not access pedestal data! No dead channel removal applied");
1327 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1330 AliFatal("Did not get geometry from EMCALLoader");
1331 return kTRUE; //coverity
1341 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1343 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1344 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1346 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1348 if (channelStatus == AliCaloCalibPedestal::kDead)