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)
781 if (!fInit) { // to prevent overwrite existing file
782 Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
786 if (strstr(option,"print")) {
792 if(strstr(option,"tim"))
793 gBenchmark->Start("EMCALDigitizer");
795 AliRunLoader *rl = AliRunLoader::Instance();
796 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
799 AliFatal("Did not get the Loader");
803 if (fLastEvent == -1) {
804 fLastEvent = rl->GetNumberOfEvents() - 1 ;
807 fLastEvent = fFirstEvent ; // what is this ??
809 nEvents = fLastEvent - fFirstEvent + 1;
812 AliEMCALGeometry *geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
813 const Int_t nTRU = geom->GetNTotalTRU();
814 TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit", nTRU*96);
815 TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
817 rl->LoadSDigits("EMCAL");
818 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
820 rl->GetEvent(ievent);
822 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
827 //-------------------------------------
830 Digits2FastOR(digitsTMP, digitsTRG);
832 WriteDigits(digitsTRG);
834 (emcalLoader->TreeD())->Fill();
836 emcalLoader->WriteDigits( "OVERWRITE");
840 digitsTRG ->Delete();
841 digitsTMP ->Delete();
843 //-------------------------------------
845 if(strstr(option,"deb"))
847 if(strstr(option,"table")) gObjectTable->Print();
849 //increment the total number of Digits per run
850 fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
855 if(strstr(option,"tim")){
856 gBenchmark->Stop("EMCALDigitizer");
857 Float_t cputime = gBenchmark->GetCpuTime("EMCALDigitizer");
858 Float_t avcputime = cputime;
859 if(nEvents==0) avcputime = 0 ;
860 AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
864 //__________________________________________________________________
865 Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
867 // Assign a smeared time to the digit, from observed distribution in LED system (?).
872 res = TMath::Sqrt(fTimeResolutionPar0 +
873 fTimeResolutionPar1/(energy*energy) );
875 // parametrization above is for ns. Convert to seconds:
879 //____________________________________________________________________________
880 void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
882 // FEE digits afterburner to produce TRG digits
883 // we are only interested in the FEE digit deposited energy
884 // to be converted later into a voltage value
886 // push the FEE digit to its associated FastOR (numbered from 0:95)
887 // TRU is in charge of summing module digits
889 AliRunLoader *runLoader = AliRunLoader::Instance();
891 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
892 if (!emcalLoader) AliFatal("Did not get the Loader");
894 const AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
896 // build FOR from simulated digits
897 // and xfer to the corresponding TRU input (mapping)
899 TClonesArray* sdigits = emcalLoader->SDigits();
901 TClonesArray *digits = (TClonesArray*)sdigits->Clone();
903 AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
905 TIter NextDigit(digits);
906 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
908 if (IsDead(digit)) continue;
910 DecalibrateTrigger(digit);
912 Int_t id = digit->GetId();
915 if (geom && geom->GetFastORIndexFromCellIndex(id , trgid))
917 AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
919 AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
923 new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
924 d = (AliEMCALDigit*)digitsTMP->At(trgid);
934 if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
936 Int_t nSamples = geom->GetNTotalTRU();
937 Int_t *timeSamples = new Int_t[nSamples];
939 NextDigit = TIter(digitsTMP);
940 while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
944 Int_t id = digit->GetId();
945 Float_t time = 50.e-9;
947 Double_t depositedEnergy = 0.;
948 for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
950 if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
952 // FIXME: Check digit time!
953 if (depositedEnergy) {
954 depositedEnergy += gRandom->Gaus(0., .08);
955 DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
957 for (Int_t j=0;j<nSamples;j++) {
958 if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
959 timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
962 new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
964 if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
969 delete [] timeSamples;
971 if (digits && digits->GetEntriesFast()) digits->Delete();
974 //____________________________________________________________________________
975 void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
979 const Int_t reso = 12; // 11-bit resolution ADC
980 const Double_t vFSR = 2.; // Full scale input voltage range 2V (p-p)
981 //const Double_t dNe = 125; // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
982 const Double_t dNe = 125/1.3; // F-ALTRO max V. FEE: factor 4
983 const Double_t vA = .136e-6; // CSP output range: 0.136uV/e-
984 const Double_t rise = 50e-9; // rise time (10-90%) of the FastOR signal before shaping
986 const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
988 Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
990 TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
991 signalF.SetParameter( 0, vV );
992 signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
993 signalF.SetParameter( 2, rise );
995 for (Int_t iTime=0; iTime<nSamples; iTime++)
997 // FIXME: add noise (probably not simply Gaussian) according to DA measurements
998 // probably plan an access to OCDB
999 Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
1000 if (TMath::Abs(sig) > vFSR/2.) {
1001 AliError("Signal overflow!");
1002 timeSamples[iTime] = (1 << reso) - 1;
1004 AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
1005 timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
1010 //____________________________________________________________________________
1011 Bool_t AliEMCALDigitizer::Init()
1013 // Makes all memory allocations
1015 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1017 if ( emcalLoader == 0 ) {
1018 Fatal("Init", "Could not obtain the AliEMCALLoader");
1023 fLastEvent = fFirstEvent ;
1026 fInput = fDigInput->GetNinputs() ;
1030 fInputFileNames = new TString[fInput] ;
1031 fEventNames = new TString[fInput] ;
1032 fInputFileNames[0] = GetTitle() ;
1033 fEventNames[0] = fEventFolderName.Data() ;
1035 for (index = 1 ; index < fInput ; index++) {
1036 fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1037 TString tempo = fDigInput->GetInputFolderName(index) ;
1038 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
1041 //Calibration instance
1042 fCalibData = emcalLoader->CalibData();
1046 //____________________________________________________________________________
1047 void AliEMCALDigitizer::InitParameters()
1049 // Parameter initialization for digitizer
1051 // Get the parameters from the OCDB via the loader
1052 AliRunLoader *rl = AliRunLoader::Instance();
1053 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1054 AliEMCALSimParam * simParam = 0x0;
1055 if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1058 simParam = AliEMCALSimParam::GetInstance();
1059 AliWarning("Simulation Parameters not available in OCDB?");
1062 fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400; // electrons per GeV
1063 fGainFluctuations = simParam->GetGainFluctuations() ; // 15.0;
1065 fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1066 if (fPinNoise < 0.0001 )
1067 Warning("InitParameters", "No noise added\n") ;
1068 fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
1069 fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1070 fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1071 fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1072 fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1074 // These defaults are normally not used.
1075 // Values are read from calibration database instead
1076 fADCchannelEC = 0.0162; // Nominal value set online for most of the channels, MIP peak sitting at ~266./16/1024
1077 fADCpedestalEC = 0.0 ; // GeV
1078 fADCchannelECDecal = 1.0; // No decalibration by default
1079 fTimeChannel = 0.0; // No time calibration by default
1080 fTimeChannelDecal = 0.0; // No time decalibration by default
1082 fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
1084 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",
1085 fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1089 //__________________________________________________________________
1090 void AliEMCALDigitizer::Print1(Option_t * option)
1091 { // 19-nov-04 - just for convenience
1093 PrintDigits(option);
1096 //__________________________________________________________________
1097 void AliEMCALDigitizer::Print (Option_t * ) const
1099 // Print Digitizer's parameters
1100 printf("Print: \n------------------- %s -------------", GetName() ) ;
1101 if( strcmp(fEventFolderName.Data(), "") != 0 ){
1102 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1106 nStreams = GetNInputStreams() ;
1113 for (index = 0 ; index < nStreams ; index++) {
1114 TString tempo(fEventNames[index]) ;
1116 if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1117 rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1118 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1120 TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1121 if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1122 fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
1123 printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1127 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1129 if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1130 else printf("\nNULL LOADER");
1132 printf("\nWith following parameters:\n") ;
1133 printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1134 printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
1135 printf("---------------------------------------------------\n") ;
1138 printf("Print: AliEMCALDigitizer not initialized") ;
1141 //__________________________________________________________________
1142 void AliEMCALDigitizer::PrintDigits(Option_t * option)
1144 //utility method for printing digit information
1146 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1148 TClonesArray * digits = emcalLoader->Digits() ;
1149 TClonesArray * sdigits = emcalLoader->SDigits() ;
1151 printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1152 printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1154 if(strstr(option,"all")){
1156 AliEMCALDigit * digit;
1157 printf("\nEMC digits (with primaries):\n") ;
1158 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1160 for (index = 0 ; index < digits->GetEntries() ; index++) {
1161 digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1163 printf("\n%6d %8f %6.5e %4d %2d : ",
1164 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1166 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1167 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1174 else printf("NULL LOADER, cannot print\n");
1177 //__________________________________________________________________
1178 Float_t AliEMCALDigitizer::TimeOfNoise(void)
1180 // Calculates the time signal generated by noise
1181 //printf("Time noise %e\n",fTimeNoise);
1182 return gRandom->Rndm() * fTimeNoise;
1185 //__________________________________________________________________
1186 void AliEMCALDigitizer::Unload()
1188 // Unloads the SDigits and Digits
1192 for(i = 1 ; i < fInput ; i++){
1193 TString tempo(fEventNames[i]) ;
1195 if ((rl = AliRunLoader::GetRunLoader(tempo)))
1196 rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1198 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1199 if(emcalLoader)emcalLoader->UnloadDigits() ;
1200 else AliFatal("Did not get the loader");
1203 //_________________________________________________________________________________________
1204 void AliEMCALDigitizer::WriteDigits()
1205 { // Makes TreeD in the output file.
1206 // Check if branch already exists:
1207 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1208 // already existing branches.
1209 // else creates branch with Digits, named "EMCAL", title "...",
1210 // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1211 // and names of files, from which digits are made.
1213 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1216 const TClonesArray * digits = emcalLoader->Digits() ;
1217 TTree * treeD = emcalLoader->TreeD();
1219 emcalLoader->MakeDigitsContainer();
1220 treeD = emcalLoader->TreeD();
1223 // -- create Digits branch
1224 Int_t bufferSize = 32000 ;
1225 TBranch * digitsBranch = 0;
1226 if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1227 digitsBranch->SetAddress(&digits);
1228 AliWarning("Digits Branch already exists. Not all digits will be visible");
1231 treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1232 //digitsBranch->SetTitle(fEventFolderName);
1236 emcalLoader->WriteDigits("OVERWRITE");
1237 emcalLoader->WriteDigitizer("OVERWRITE");
1243 else AliFatal("Loader not available");
1246 //__________________________________________________________________
1247 void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1248 { // overloaded method
1249 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1252 TTree* treeD = emcalLoader->TreeD();
1255 emcalLoader->MakeDigitsContainer();
1256 treeD = emcalLoader->TreeD();
1259 // -- create Digits branch
1260 Int_t bufferSize = 32000;
1262 if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1264 triggerBranch->SetAddress(&digits);
1268 treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1273 else AliFatal("Loader not available");
1276 //__________________________________________________________________
1277 Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1279 AliRunLoader *runLoader = AliRunLoader::Instance();
1280 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1281 if (!emcalLoader) AliFatal("Did not get the Loader");
1283 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1285 AliWarning("Could not access pedestal data! No dead channel removal applied");
1290 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1291 if (!geom) AliFatal("Did not get geometry from EMCALLoader");
1293 Int_t absId = digit->GetId();
1301 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1303 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1304 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1306 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1308 if (channelStatus == AliCaloCalibPedestal::kDead)
1315 //__________________________________________________________________
1316 Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1318 AliRunLoader *runLoader = AliRunLoader::Instance();
1319 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1320 if (!emcalLoader) AliFatal("Did not get the Loader");
1322 AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1324 AliWarning("Could not access pedestal data! No dead channel removal applied");
1329 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1330 if (!geom) AliFatal("Did not get geometry from EMCALLoader");
1339 Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1341 if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1342 geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1344 Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1346 if (channelStatus == AliCaloCalibPedestal::kDead)