* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-
/* $Id$ */
-
+/** @file AliFMDDigitizer.cxx
+ @author Christian Holm Christensen <cholm@nbi.dk>
+ @date Mon Mar 27 12:38:26 2006
+ @brief FMD Digitizers implementation
+ @ingroup FMD_sim
+*/
//////////////////////////////////////////////////////////////////////////////
//
// This class contains the procedures simulation ADC signal for the
//
#include <TTree.h> // ROOT_TTree
-#include <TRandom.h> // ROOT_TRandom
+//#include <TRandom.h> // ROOT_TRandom
#include <AliLog.h> // ALILOG_H
#include "AliFMDDigitizer.h" // ALIFMDDIGITIZER_H
#include "AliFMD.h" // ALIFMD_H
-#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
-#include "AliFMDDetector.h" // ALIFMDDETECTOR_H
-#include "AliFMDRing.h" // ALIFMDRING_H
-#include "AliFMDHit.h" // ALIFMDHIT_H
+// #include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
+// #include "AliFMDDetector.h" // ALIFMDDETECTOR_H
+// #include "AliFMDRing.h" // ALIFMDRING_H
+// #include "AliFMDHit.h" // ALIFMDHIT_H
#include "AliFMDDigit.h" // ALIFMDDIGIT_H
+#include "AliFMDParameters.h" // ALIFMDPARAMETERS_H
#include <AliRunDigitizer.h> // ALIRUNDIGITIZER_H
#include <AliRun.h> // ALIRUN_H
#include <AliLoader.h> // ALILOADER_H
#include <AliRunLoader.h> // ALIRUNLOADER_H
-//====================================================================
-ClassImp(AliFMDBaseDigitizer)
-#if 0
- ; // This is here to keep Emacs for indenting the next line
-#endif
-
-//____________________________________________________________________
-AliFMDBaseDigitizer::AliFMDBaseDigitizer()
- : fRunLoader(0)
-{
- // Default ctor - don't use it
-}
-
-//____________________________________________________________________
-AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager)
- : AliDigitizer(manager, "AliFMDBaseDigitizer", "FMD Digitizer base class"),
- fRunLoader(0),
- fEdep(AliFMDMap::kMaxDetectors,
- AliFMDMap::kMaxRings,
- AliFMDMap::kMaxSectors,
- AliFMDMap::kMaxStrips)
-{
- // Normal CTOR
- AliDebug(1," processed");
- SetVA1MipRange();
- SetAltroChannelSize();
- SetSampleRate();
- SetShapingTime();
-}
-
-//____________________________________________________________________
-AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name,
- const Char_t* title)
- : AliDigitizer(name, title),
- fRunLoader(0),
- fEdep(AliFMDMap::kMaxDetectors,
- AliFMDMap::kMaxRings,
- AliFMDMap::kMaxSectors,
- AliFMDMap::kMaxStrips)
-{
- // Normal CTOR
- AliDebug(1," processed");
- SetVA1MipRange();
- SetAltroChannelSize();
- SetSampleRate();
- SetShapingTime();
-}
-
-//____________________________________________________________________
-AliFMDBaseDigitizer::~AliFMDBaseDigitizer()
-{
- // Destructor
-}
-
-//____________________________________________________________________
-Bool_t
-AliFMDBaseDigitizer::Init()
-{
- // Initialization
- return kTRUE;
-}
-
-
-//____________________________________________________________________
-void
-AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
-{
- // Sum energy deposited contributions from each hit in a cache
- // (fEdep).
- if (!fRunLoader)
- Fatal("SumContributions", "no run loader");
-
- // Clear array of deposited energies
- fEdep.Reset();
-
- // Get the FMD loader
- AliLoader* inFMD = fRunLoader->GetLoader("FMDLoader");
- // And load the hits
- inFMD->LoadHits("READ");
-
- // Get the tree of hits
- TTree* hitsTree = inFMD->TreeH();
- if (!hitsTree) {
- // Try again
- inFMD->LoadHits("READ");
- hitsTree = inFMD->TreeH();
- }
-
- // Get the FMD branch
- TBranch* hitsBranch = hitsTree->GetBranch("FMD");
- if (hitsBranch) fmd->SetHitsAddressBranch(hitsBranch);
- else AliFatal("Branch FMD hit not found");
-
- // Get a list of hits from the FMD manager
- TClonesArray *fmdHits = fmd->Hits();
-
- // Get number of entries in the tree
- Int_t ntracks = Int_t(hitsTree->GetEntries());
-
- // Loop over the tracks in the
- for (Int_t track = 0; track < ntracks; track++) {
- // Read in entry number `track'
- hitsBranch->GetEntry(track);
-
- // Get the number of hits
- Int_t nhits = fmdHits->GetEntries ();
- for (Int_t hit = 0; hit < nhits; hit++) {
- // Get the hit number `hit'
- AliFMDHit* fmdHit =
- static_cast<AliFMDHit*>(fmdHits->UncheckedAt(hit));
-
- // Extract parameters
- UShort_t detector = fmdHit->Detector();
- Char_t ring = fmdHit->Ring();
- UShort_t sector = fmdHit->Sector();
- UShort_t strip = fmdHit->Strip();
- Float_t edep = fmdHit->Edep();
- if (fEdep(detector, ring, sector, strip).fEdep != 0)
- AliDebug(1, Form("Double hit in %d%c(%d,%d)",
- detector, ring, sector, strip));
-
- fEdep(detector, ring, sector, strip).fEdep += edep;
- fEdep(detector, ring, sector, strip).fN += 1;
- // Add this to the energy deposited for this strip
- } // hit loop
- } // track loop
-}
-
-//____________________________________________________________________
-void
-AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
-{
- // For the stored energy contributions in the cache (fEdep), convert
- // the energy signal to ADC counts, and store the created digit in
- // the digits array (AliFMD::fDigits)
- //
- AliFMDGeometry* geometry = AliFMDGeometry::Instance();
-
- TArrayI counts(3);
- for (UShort_t detector=1; detector <= 3; detector++) {
- // Get pointer to subdetector
- AliFMDDetector* det = geometry->GetDetector(detector);
- if (!det) continue;
- for (UShort_t ringi = 0; ringi <= 1; ringi++) {
- // Get pointer to Ring
- AliFMDRing* r = det->GetRing((ringi == 0 ? 'I' : 'O'));
- if (!r) continue;
-
- // Get number of sectors
- UShort_t nSectors = UShort_t(360. / r->GetTheta());
- // Loop over the number of sectors
- for (UShort_t sector = 0; sector < nSectors; sector++) {
- // Get number of strips
- UShort_t nStrips = r->GetNStrips();
- // Loop over the stips
- Float_t last = 0;
- for (UShort_t strip = 0; strip < nStrips; strip++) {
- // Reset the counter array to the invalid value -1
- counts.Reset(-1);
- // Reset the last `ADC' value when we've get to the end of a
- // VA1_ALICE channel.
- if (strip % 128 == 0) last = 0;
-
- Float_t edep = fEdep(detector, r->GetId(), sector, strip).fEdep;
- ConvertToCount(edep, last, r->GetSiThickness(),
- geometry->GetSiDensity(), counts);
- last = edep;
- AddDigit(fmd, detector, r->GetId(), sector, strip,
- edep, UShort_t(counts[0]),
- Short_t(counts[1]), Short_t(counts[2]));
-#if 0
- // This checks if the digit created will give the `right'
- // number of particles when reconstructed, using a naiive
- // approach. It's here only as a quality check - nothing
- // else.
- CheckDigit(fEdep(detector, r->GetId(), sector, strip).fEdep,
- fEdep(detector, r->GetId(), sector, strip).fN,
- counts);
-#endif
- } // Strip
- } // Sector
- } // Ring
- } // Detector
-}
-
-//____________________________________________________________________
-void
-AliFMDBaseDigitizer::ConvertToCount(Float_t edep,
- Float_t last,
- Float_t siThickness,
- Float_t siDensity,
- TArrayI& counts) const
-{
- // Convert the total energy deposited to a (set of) ADC count(s).
- //
- // This is done by
- //
- // Energy_Deposited ALTRO_Channel_Size
- // ADC = -------------------------- ------------------- + pedestal
- // Energy_Deposition_Of_1_MIP VA1_ALICE_MIP_Range
- //
- // Energy_Deposited fAltroChannelSize
- // = --------------------------------- ----------------- + pedestal
- // 1.664 * Si_Thickness * Si_Density fVA1MipRange
- //
- //
- // = Energy_Deposited * ConversionFactor + pedestal
- //
- // However, this is modified by the response function of the
- // VA1_ALICE pre-amp. chip in case we are doing oversampling of the
- // VA1_ALICE output.
- //
- // In that case, we get N=fSampleRate values of the ADC, and the
- // `EnergyDeposited' is a function of which sample where are
- // calculating the ADC for
- //
- // ADC_i = f(EnergyDeposited, i/N, Last) * ConversionFactor + pedestal
- //
- // where Last is the Energy deposited in the previous strip.
- //
- // Here, f is the shaping function of the VA1_ALICE. This is given
- // by
- //
- // | (E - l) * (1 - exp(-B * t) + l if E > l
- // f(E, t, l) = <
- // | (l - E) * exp(-B * t) + E otherwise
- //
- //
- // = E + (l - E) * ext(-B * t)
- //
- const Float_t mipEnergy = 1.664 * siThickness * siDensity;
- const Float_t convf = (1 / mipEnergy * Float_t(fAltroChannelSize)
- / fVA1MipRange);
- UShort_t ped = MakePedestal();
-
- // In case we don't oversample, just return the end value.
- if (fSampleRate == 1) {
- counts[0] = UShort_t(TMath::Min(edep * convf + ped,
- Float_t(fAltroChannelSize)));
- return;
- }
-
- // Create a pedestal
- Int_t n = fSampleRate;
- Float_t B = fShapingTime;
- for (Ssiz_t i = 0; i < n; i++) {
- Float_t t = Float_t(i) / n;
- Float_t s = edep + (last - edep) * TMath::Exp(-B * t);
- counts[i] = UShort_t(TMath::Min(s * convf + ped,
- Float_t(fAltroChannelSize)));
- }
-}
-
-
//====================================================================
ClassImp(AliFMDDigitizer)
{
// Normal CTOR
AliDebug(1," processed");
- SetPedestal();
}
//____________________________________________________________________
if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
// Get the AliFMD object
- AliFMD* fmd = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
+ AliFMD* fmd =
+ static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
if (!fmd) {
AliError("Can not get FMD from gAlice");
return;
fmd->MakeBranchInTree(digitTree, fmd->GetName(), &(digits), 4000, 0);
// TBranch* digitBranch = digitTree->GetBranch(fmd->GetName());
// Fill the tree
- digitTree->Fill();
+ Int_t write = 0;
+ write = digitTree->Fill();
+ AliDebug(1, Form("Wrote %d bytes to digit tree", write));
// Write the digits to disk
outFMD->WriteDigits("OVERWRITE");
//____________________________________________________________________
UShort_t
-AliFMDDigitizer::MakePedestal() const
+AliFMDDigitizer::MakePedestal(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const
{
- return UShort_t(TMath::Max(gRandom->Gaus(fPedestal, fPedestalWidth), 0.));
+ // Make a pedestal
+ AliFMDParameters* param =AliFMDParameters::Instance();
+ Float_t mean =param->GetPedestal(detector,ring,sector,strip);
+ Float_t width =param->GetPedestalWidth(detector,ring,sector,strip);
+ return UShort_t(TMath::Max(gRandom->Gaus(mean, width), 0.));
}
//____________________________________________________________________
Short_t count2,
Short_t count3) const
{
+ // Add a digit
fmd->AddDigitByFields(detector, ring, sector, strip, count1, count2, count3);
}
//____________________________________________________________________
void
-AliFMDDigitizer::CheckDigit(Float_t /* edep */,
+AliFMDDigitizer::CheckDigit(AliFMDDigit* digit,
UShort_t nhits,
const TArrayI& counts)
{
- Int_t integral = counts[0];
+ // Check that digit is consistent
+ AliFMDParameters* param = AliFMDParameters::Instance();
+ UShort_t det = digit->Detector();
+ Char_t ring = digit->Ring();
+ UShort_t sec = digit->Sector();
+ UShort_t str = digit->Strip();
+ Float_t mean = param->GetPedestal(det,ring,sec,str);
+ Float_t width = param->GetPedestalWidth(det,ring,sec,str);
+ UShort_t range = param->GetVA1MipRange();
+ UShort_t size = param->GetAltroChannelSize();
+ Int_t integral = counts[0];
if (counts[1] >= 0) integral += counts[1];
if (counts[2] >= 0) integral += counts[2];
- integral -= Int_t(fPedestal + 2 * fPedestalWidth);
+ integral -= Int_t(mean + 2 * width);
if (integral < 0) integral = 0;
- Float_t convf = Float_t(fVA1MipRange) / fAltroChannelSize;
- Float_t mips = integral * convf;
+ Float_t convF = Float_t(range) / size;
+ Float_t mips = integral * convF;
if (mips > Float_t(nhits) + .5 || mips < Float_t(nhits) - .5)
Warning("CheckDigit", "Digit -> %4.2f MIPS != %d +/- .5 hits",
mips, nhits);
}
-//====================================================================
-ClassImp(AliFMDSDigitizer)
-
-//____________________________________________________________________
-AliFMDSDigitizer::AliFMDSDigitizer()
-{
- // Default ctor - don't use it
-}
-
-//____________________________________________________________________
-AliFMDSDigitizer::AliFMDSDigitizer(const Char_t* headerFile,
- const Char_t* /* sdigfile */)
- : AliFMDBaseDigitizer("FMDSDigitizer", "FMD SDigitizer")
-{
- // Normal CTOR
- AliDebug(1," processed");
-
- fRunLoader = AliRunLoader::GetRunLoader(); // Open(headerFile);
- if (!fRunLoader)
- Fatal("AliFMDSDigitizer", "cannot open session, header file '%s'",
- headerFile);
- AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
- if (!loader)
- Fatal("AliFMDSDigitizer", "cannot find FMD loader in specified event");
-
- // Add task to tasks folder
- loader->PostSDigitizer(this);
-
-}
-
-//____________________________________________________________________
-AliFMDSDigitizer::~AliFMDSDigitizer()
-{
- AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
- loader->CleanSDigitizer();
-}
-
-//____________________________________________________________________
-void
-AliFMDSDigitizer::Exec(Option_t*)
-{
- // Get the output manager
- if (!fRunLoader) {
- Error("Exec", "Run loader is not set");
- return;
- }
- if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
- if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
-
- AliLoader* fmdLoader = fRunLoader->GetLoader("FMDLoader");
- if (!fmdLoader) Fatal("Exec", "no FMD loader");
-
- // Get the AliFMD object
- AliFMD* fmd =
- static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
- if (!fmd) {
- AliError("Can not get FMD from gAlice");
- return;
- }
-
- Int_t nEvents = Int_t(fRunLoader->TreeE()->GetEntries());
- for (Int_t event = 0; event < nEvents; event++) {
- AliDebug(1,Form(" Digitizing event number %d", event));
- // Get the current loader
- fRunLoader->GetEvent(event);
-
- if (!fmdLoader->TreeS()) fmdLoader->MakeTree("S");
- // Make a branch
- fmd->MakeBranch("S");
-
- // Cache contriutions
- SumContributions(fmd);
-
- // Digitize the event
- DigitizeHits(fmd);
-
- fmdLoader->TreeS()->Reset();
- fmdLoader->TreeS()->Fill();
- fmdLoader->WriteSDigits("OVERWRITE");
- }
-}
-
-//____________________________________________________________________
-void
-AliFMDSDigitizer::AddDigit(AliFMD* fmd,
- UShort_t detector,
- Char_t ring,
- UShort_t sector,
- UShort_t strip,
- Float_t edep,
- UShort_t count1,
- Short_t count2,
- Short_t count3) const
-{
- fmd->AddSDigitByFields(detector, ring, sector, strip, edep,
- count1, count2, count3);
-}
-
-
-
//____________________________________________________________________
//
// EOF