/**************************************************************************
- * Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
+ * Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
* *
* Author: The ALICE Off-line Project. *
* Contributors are mentioned in the code where appropriate. *
* 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
+// Forward Multiplicity detector : SDigits->Digits
+//
+// Digits consists of
+// - Detector #
+// - Ring ID
+// - Sector #
+// - Strip #
+// - ADC count in this channel
+//
+// Digits consists of
+// - Detector #
+// - Ring ID
+// - Sector #
+// - Strip #
+// - Total energy deposited in the strip
+// - ADC count in this channel
+//
+// As the Digits and SDigits have so much in common, the classes
+// AliFMDDigitizer and AliFMDSDigitizer are implemented via a base
+// class AliFMDBaseDigitizer.
+//
+// +---------------------+
+// | AliFMDBaseDigitizer |
+// +---------------------+
+// ^
+// |
+// +----------+---------+
+// | |
+// +-----------------+ +------------------+
+// | AliFMDDigitizer | | AliFMDSDigitizer |
+// +-----------------+ +------------------+
+// |
+// +-------------------+
+// | AliFMDSSDigitizer |
+// +-------------------+
+//
+// These classes has several paramters:
+//
+// fPedestal
+// fPedestalWidth
+// (Only AliFMDDigitizer)
+// Mean and width of the pedestal. The pedestal is simulated
+// by a Guassian, but derived classes my override MakePedestal
+// to simulate it differently (or pick it up from a database).
+//
+// fVA1MipRange
+// The dymamic MIP range of the VA1_ALICE pre-amplifier chip
+//
+// fAltroChannelSize
+// The largest number plus one that can be stored in one
+// channel in one time step in the ALTRO ADC chip.
+//
+// fSampleRate
+// How many times the ALTRO ADC chip samples the VA1_ALICE
+// pre-amplifier signal. The VA1_ALICE chip is read-out at
+// 10MHz, while it's possible to drive the ALTRO chip at
+// 25MHz. That means, that the ALTRO chip can have time to
+// sample each VA1_ALICE signal up to 2 times. Although it's
+// not certain this feature will be used in the production,
+// we'd like have the option, and so it should be reflected in
+// the code.
+//
+//
+// The shaping function of the VA1_ALICE is generally given by
+//
+// f(x) = A(1 - exp(-Bx))
+//
+// where A is the total charge collected in the pre-amp., and B is a
+// paramter that depends on the shaping time of the VA1_ALICE circut.
+//
+// When simulating the shaping function of the VA1_ALICe
+// pre-amp. chip, we have to take into account, that the shaping
+// function depends on the previous value of read from the pre-amp.
+//
+// That results in the following algorithm:
+//
+// last = 0;
+// FOR charge IN pre-amp. charge train DO
+// IF last < charge THEN
+// f(t) = (charge - last) * (1 - exp(-B * t)) + last
+// ELSE
+// f(t) = (last - charge) * exp(-B * t) + charge)
+// ENDIF
+// FOR i IN # samples DO
+// adc_i = f(i / (# samples))
+// DONE
+// last = charge
+// DONE
+//
+// Here,
+//
+// pre-amp. charge train
+// is a series of 128 charges read from the VA1_ALICE chip
+//
+// # samples
+// is the number of times the ALTRO ADC samples each of the 128
+// charges from the pre-amp.
+//
+// Where Q is the total charge collected by the VA1_ALICE
+// pre-amplifier. Q is then given by
+//
+// E S
+// Q = - -
+// e R
+//
+// where E is the total energy deposited in a silicon strip, R is the
+// dynamic range of the VA1_ALICE pre-amp (fVA1MipRange), e is the
+// energy deposited by a single MIP, and S ALTRO channel size in each
+// time step (fAltroChannelSize).
+//
+// The energy deposited per MIP is given by
+//
+// e = M * rho * w
+//
+// where M is the universal number 1.664, rho is the density of
+// silicon, and w is the depth of the silicon sensor.
+//
+// The final ADC count is given by
+//
+// C' = C + P
+//
+// where P is the (randomized) pedestal (see MakePedestal)
+//
+// This class uses the class template AliFMDMap<Type> to make an
+// internal cache of the energy deposted of the hits. The class
+// template is instantasized as
+//
+// typedef AliFMDMap<std::pair<Float_t, UShort_t> > AliFMDEdepMap;
+//
+// The first member of the values is the summed energy deposition in a
+// given strip, while the second member of the values is the number of
+// hits in a given strip. Using the second member, it's possible to
+// do some checks on just how many times a strip got hit, and what
+// kind of error we get in our reconstructed hits. Note, that this
+// information is currently not written to the digits tree. I think a
+// QA (Quality Assurance) digit tree is better suited for that task.
+// However, the information is there to be used in the future.
+//
+//
+// Latest changes by Christian Holm Christensen
+//
+//////////////////////////////////////////////////////////////////////////////
+// /1
+// | A(-1 + B + exp(-B))
+// | f(x) dx = ------------------- = 1
+// | B
+// / 0
+//
+// and B is the a parameter defined by the shaping time (fShapingTime).
+//
+// Solving the above equation, for A gives
+//
+// B
+// A = ----------------
+// -1 + B + exp(-B)
+//
+// So, if we define the function g: [0,1] -> [0:1] by
+//
+// / v
+// | Bu + exp(-Bu) - Bv - exp(-Bv)
+// g(u,v) = | f(x) dx = -A -----------------------------
+// | B
+// / u
+//
+// we can evaluate the ALTRO sample of the VA1_ALICE pre-amp between
+// any two times (u, v), by
+//
+//
+// B Bu + exp(-Bu) - Bv - exp(-Bv)
+// C = Q g(u,v) = - Q ---------------- -----------------------------
+// -1 + B + exp(-B) B
+//
+// Bu + exp(-Bu) - Bv - exp(-Bv)
+// = - Q -----------------------------
+// -1 + B + exp(-B)
+//
-#include <TTree.h>
-#include <TVector.h>
-#include <TObjArray.h>
+#include <TTree.h> // ROOT_TTree
#include <TFile.h>
-#include <TDirectory.h>
-#include <TRandom.h>
+#include "AliFMDDebug.h" // Better debug macros
+#include "AliFMDDigitizer.h" // ALIFMDSSDIGITIZER_H
+#include "AliFMD.h" // ALIFMD_H
+#include "AliFMDSDigit.h" // ALIFMDDIGIT_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(AliFMDDigitizer)
+#if 0
+;
+#endif
+//____________________________________________________________________
+Bool_t
+AliFMDDigitizer::Init()
+{
+ // Initialisation
+ if (!AliFMDBaseDigitizer::Init()) return kFALSE;
+
+#if 0
+ // Get the AliRun object
+ AliRun* run = fRunLoader->GetAliRun();
+ if (!run) {
+ AliWarning("Loading gAlice");
+ fRunLoader->LoadgAlice();
+ if (!run) {
+ AliError("Can not get Run from Run Loader");
+ return kFALSE;
+ }
+ }
+
+ // Get the AliFMD object
+ fFMD = static_cast<AliFMD*>(run->GetDetector("FMD"));
+ if (!fFMD) {
+ AliError("Can not get FMD from gAlice");
+ return kFALSE;
+ }
+#endif
+ return kTRUE;
+}
-#include "AliFMDDigitizer.h"
-#include "AliFMD.h"
-#include "AliFMDSDigitizer.h"
-#include "AliFMDhit.h"
-#include "AliFMDdigit.h"
-#include "AliRunDigitizer.h"
-#include "AliRun.h"
-#include "AliPDG.h"
+//____________________________________________________________________
+void
+AliFMDDigitizer::Exec(Option_t*)
+{
+ if (!fManager) {
+ AliError("No digitisation manager defined");
+ return;
+ }
-#include <stdlib.h>
-#include <Riostream.h>
-#include <Riostream.h>
+ // Clear array of deposited energies
+ fEdep.Reset();
-ClassImp(AliFMDDigitizer)
+ AliRunLoader* runLoader = 0;
+ if (!gAlice) {
+ TString folderName(fManager->GetInputFolderName(0));
+ runLoader = AliRunLoader::GetRunLoader(folderName.Data());
+ if (!runLoader) {
+ AliError(Form("Failed at getting run loader from %s",
+ folderName.Data()));
+ return;
+ }
+ if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
+ runLoader->GetAliRun();
+ }
+ if (!gAlice) {
+ AliError("Can not get Run from Run Loader");
+ return;
+ }
+
+ // Get the AliFMD object
+ fFMD = static_cast<AliFMD*>(gAlice->GetDetector("FMD"));
+ if (!fFMD) {
+ AliError("Can not get FMD from gAlice");
+ return;
+ }
-//___________________________________________
- AliFMDDigitizer::AliFMDDigitizer() :AliDigitizer()
-{
-// Default ctor - don't use it
- ;
-}
-//___________________________________________
-AliFMDDigitizer::AliFMDDigitizer(AliRunDigitizer* manager)
- :AliDigitizer(manager)
-{
- cout<<"AliFMDDigitizer::AliFMDDigitizer"<<endl;
-// ctor which should be used
-// fDebug =0;
- // if (GetDebug()>2)
- // cerr<<"AliFMDDigitizer::AliFMDDigitizer"
- // <<"(AliRunDigitizer* manager) was processed"<<endl;
-}
+ // Loop over input files
+ Int_t nFiles= fManager->GetNinputs();
+ AliFMDDebug(1, (" Digitizing event number %d, got %d inputs",
+ fManager->GetOutputEventNr(), nFiles));
+ for (Int_t inputFile = 0; inputFile < nFiles; inputFile++) {
+ AliFMDDebug(5, ("Now reading input # %d", inputFile));
+ // Get the current loader
+ AliRunLoader* currentLoader =
+ AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
+ if (!currentLoader) {
+ Error("Exec", Form("no run loader for input file # %d", inputFile));
+ continue;
+ }
-//------------------------------------------------------------------------
-AliFMDDigitizer::~AliFMDDigitizer()
-{
-// Destructor
-}
+ // Cache contriutions
+ AliFMDDebug(5, ("Now summing the contributions from input # %d",inputFile));
- //------------------------------------------------------------------------
-/*
-Bool_t AliFMDDigitizer::Init()
-{
- cout<<"AliFMDDigitizer::Init"<<endl;
- return kTRUE;
-}
-
-*/
-//---------------------------------------------------------------------
+ // Get the FMD loader
+ AliLoader* inFMD = currentLoader->GetLoader("FMDLoader");
+ // And load the summable digits
+ inFMD->LoadSDigits("READ");
+
+ // Get the tree of summable digits
+ TTree* sdigitsTree = inFMD->TreeS();
+ if (!sdigitsTree) {
+ AliError("No sdigit tree from manager");
+ continue;
+ }
+ if (AliLog::GetDebugLevel("FMD","") >= 10) {
+ TFile* file = sdigitsTree->GetCurrentFile();
+ if (!file) {
+ AliWarning("Input tree has no file!");
+ }
+ else {
+ AliFMDDebug(10, ("Input tree file %s content:", file->GetName()));
+ file->ls();
+ }
+ // AliFMDDebug(5, ("Input tree %s file structure:",
+ // sdigitsTree->GetName()));
+ // sdigitsTree->Print();
+ }
-void AliFMDDigitizer::Exec(Option_t* option)
-{
- /*
- Conver hits to digits:
- - number of detector;
- - number of ring;
- - number of sector;
- - ADC signal in this channel
- */
+ // Get the FMD branch
+ TBranch* sdigitsBranch = sdigitsTree->GetBranch("FMD");
+ if (!sdigitsBranch) {
+ AliError("Failed to get sdigit branch");
+ return;
+ }
+ // Set the branch addresses
+ fFMD->SetTreeAddress();
-#ifdef DEBUG
- cout<<"AliFMDDigitizer::>SDigits2Digits start...\n";
-#endif
+ // Sum contributions from the sdigits
+ AliFMDDebug(3, ("Will now sum contributions from SDigits"));
+ SumContributions(sdigitsBranch);
+
+ // Unload the sdigits
+ inFMD->UnloadSDigits();
+ }
- // cout<<" FMD "<<FMD<<endl;
-
- Int_t volume, sector, ring, charge;
- Float_t e;
- Float_t de[10][50][800];
- Int_t hit;
- Int_t digit[5];
- Int_t ivol, iSector, iRing;
- for (Int_t i=0; i<10; i++)
- for(Int_t j=0; j<50; j++)
- for(Int_t ij=0; ij<800; ij++)
- de[i][j][ij]=0;
- Int_t numberOfRings[5]=
- {768,384,768,384,768};
- Int_t numberOfSector[5]=
- {20,40,20,40,20};
+ TString outFolder(fManager->GetOutputFolderName());
+ AliRunLoader* out = AliRunLoader::GetRunLoader(outFolder.Data());
+ AliLoader* outFMD = out->GetLoader("FMDLoader");
+ if (!outFMD) {
+ AliError("Cannot get the FMDLoader output folder");
+ return;
+ }
+ TTree* outTree = MakeOutputTree(outFMD);
+ if (!outTree) {
+ AliError("Failed to get output tree");
+ return;
+ }
+ // Set the branch address
+ fFMD->SetTreeAddress();
- AliFMDhit *fmdHit=0;
- TTree *tH=0;
- TBranch *brHits=0;
-
- AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD") ;
-
-// Loop over files to digitize
-
- Int_t nFiles=GetManager()->GetNinputs();
- for (Int_t inputFile=0; inputFile<nFiles;
- inputFile++) {
-
- cout<<" event "<<fManager->GetOutputEventNr()<<endl;
- if (fFMD)
- {
- TClonesArray *FMDhits = fFMD->Hits ();
- tH = fManager->GetInputTreeH(inputFile);
- brHits = tH->GetBranch("FMD");
- if (brHits) {
- fFMD->SetHitsAddressBranch(brHits);
- }else{
- cerr<<"EXEC Branch FMD hit not found"<<endl;
- exit(111);
- }
- Int_t ntracks = (Int_t) tH->GetEntries();
-
- for (Int_t track = 0; track < ntracks; track++)
- {
- brHits->GetEntry(track);
- Int_t nhits = FMDhits->GetEntries ();
-
- for (hit = 0; hit < nhits; hit++)
- {
- fmdHit = (AliFMDhit *) FMDhits->UncheckedAt (hit);
-
- volume = fmdHit->Volume ();
- sector = fmdHit->NumberOfSector ();
- ring = fmdHit->NumberOfRing ();
- e = fmdHit->Edep ();
- de[volume][sector][ring] += e;
-
- } //hit loop
- } //track loop
- } //if FMD
-
-
- // Put noise and make ADC signal
- Float_t iP = 1.664 * 0.04 * 2.33 / 22400; // = 0.69e-6;
- for ( ivol=1; ivol<=5; ivol++){
- for ( iSector=1; iSector<=numberOfSector[ivol-1]; iSector++){
- for ( iRing=1; iRing<=numberOfRings[ivol-1]; iRing++){
- digit[0]=ivol;
- digit[1]=iSector;
- digit[2]=iRing;
- charge = Int_t (de[ivol][iSector][iRing] / iP);
- digit[3]=PutNoise(charge);
- if(digit[3]<= 500) digit[3]=500;
- //dynamic range from MIP(0.155MeV) to 30MIP(4.65MeV)
- //1024 ADC channels
- Float_t channelWidth=(22400*50)/1024;
- digit[4]=Int_t(digit[3]/channelWidth);
- if (digit[4]>1024) digit[4]=1024;
- fFMD->AddDigit(digit);
- } //ivol
- } //iSector
- } //iRing
-
- TTree* treeD = fManager->GetTreeD();
- treeD->Clear();
- treeD->Reset();
- fFMD->MakeBranchInTreeD(treeD);
- treeD->Fill();
-
- fManager->GetTreeD()->Write(0,TObject::kOverwrite);
+ // And digitize the cached data
+ DigitizeHits();
- gAlice->ResetDigits();
- }
+ // Fill the tree
+ Int_t write = outTree->Fill();
+ AliFMDDebug(5, ("Wrote %d bytes to digit tree", write));
+
+ // Store the digits
+ StoreDigits(outFMD);
+}
+
+//____________________________________________________________________
+void
+AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch)
+{
+ // Sum energy deposited contributions from each sdigits in a cache
+ // (fEdep).
+ AliFMDDebug(3, ("Runnin our version of SumContributions"));
+
+ // Get a list of hits from the FMD manager
+ TClonesArray *fmdSDigits = fFMD->SDigits();
+
+ // Get number of entries in the tree
+ Int_t nevents = Int_t(sdigitsBranch->GetEntries());
+
+ Int_t read = 0;
+ // Loop over the events in the
+ for (Int_t event = 0; event < nevents; event++) {
+ // Read in entry number `event'
+ read += sdigitsBranch->GetEntry(event);
+
+ // Get the number of sdigits
+ Int_t nsdigits = fmdSDigits->GetEntries ();
+ AliFMDDebug(3, ("Got %5d SDigits", nsdigits));
+ for (Int_t sdigit = 0; sdigit < nsdigits; sdigit++) {
+ // Get the sdigit number `sdigit'
+ AliFMDSDigit* fmdSDigit =
+ static_cast<AliFMDSDigit*>(fmdSDigits->UncheckedAt(sdigit));
+
+ AliFMDDebug(5, ("Adding contribution of %d tracks",
+ fmdSDigit->GetNTrack()));
+ AliFMDDebug(15, ("Contrib from FMD%d%c[%2d,%3d] (%s) from track %d",
+ fmdSDigit->Detector(),
+ fmdSDigit->Ring(),
+ fmdSDigit->Sector(),
+ fmdSDigit->Strip(),
+ fmdSDigit->GetName(),
+ fmdSDigit->GetTrack(0)));
+
+ // Extract parameters
+ AddContribution(fmdSDigit->Detector(),
+ fmdSDigit->Ring(),
+ fmdSDigit->Sector(),
+ fmdSDigit->Strip(),
+ fmdSDigit->Edep(),
+ kTRUE,
+ fmdSDigit->GetNTrack(),
+ fmdSDigit->GetTracks());
+ } // sdigit loop
+ } // event loop
+
+
+ AliFMDDebug(3, ("Size of cache: %d bytes, read %d bytes",
+ sizeof(fEdep), read));
}
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
+