From: schutz Date: Fri, 16 Jul 2004 17:00:04 +0000 (+0000) Subject: Made the raw data format conversion from digits more modular (+2 more methods) to... X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=046ae904fb5d34b39750ee5f05225b513238ac79 Made the raw data format conversion from digits more modular (+2 more methods) to isolate the sampled response function --- diff --git a/EMCAL/AliEMCAL.cxx b/EMCAL/AliEMCAL.cxx index 3b6dd1da797..1e95f60c6b5 100644 --- a/EMCAL/AliEMCAL.cxx +++ b/EMCAL/AliEMCAL.cxx @@ -30,6 +30,9 @@ class TFile; #include #include #include +#include +#include +#include // --- Standard library --- @@ -47,17 +50,19 @@ ClassImp(AliEMCAL) AliEMCAL::AliEMCAL():AliDetector() { // Default ctor - fName="EMCAL"; - // fGeom = 0 ; - fRan = new TRandom(123456) ; + fName = "EMCAL" ; } //____________________________________________________________________________ AliEMCAL::AliEMCAL(const char* name, const char* title): AliDetector(name,title) { // ctor : title is used to identify the layout - // fGeom = 0; - fRan = new TRandom(123456) ; + + fTimeMax = 1.28E-5 ; + fTimePeak = 2.0E-6 ; + fTimeRes = 1.5E-6 ; + fHighGainFactor = 40; + fHighGainOffset = 0x200 ; } //____________________________________________________________________________ @@ -192,20 +197,13 @@ void AliEMCAL::Digits2Raw() // some digitization constants const Int_t kDDLOffset = 0x800; - const Double_t kTimeMax = 1.28E-5; - const Int_t kTimeBins = 256; - const Double_t kTimePeak = 2.0E-6; - const Double_t kTimeRes = 1.5E-6; const Int_t kThreshold = 3; - const Int_t kHighGainFactor = 40; - const Int_t kHighGainOffset = 0x200; - // PHOS has 4 DDL per module; I assume therefore that kChannelsperDDL=896+1 EMCAL channel go to one DDL const Int_t kChannelsperDDL = 897 ; AliAltroBuffer* buffer = NULL; Int_t prevDDL = -1; - Int_t adcValuesLow[kTimeBins]; - Int_t adcValuesHigh[kTimeBins]; + Int_t adcValuesLow[fkTimeBins]; + Int_t adcValuesHigh[fkTimeBins]; // loop over digits (assume ordered digits) for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { @@ -235,9 +233,9 @@ void AliEMCAL::Digits2Raw() } // out of time range signal (?) - if (digit->GetTime() > kTimeMax) { + if (digit->GetTimeR() > fTimeMax) { buffer->FillBuffer(digit->GetAmp()); - buffer->FillBuffer(kTimeBins); // time bin + buffer->FillBuffer(fkTimeBins); // time bin buffer->FillBuffer(3); // bunch length buffer->WriteTrailer(3, idDDL, 0, 0); // trailer @@ -245,48 +243,21 @@ void AliEMCAL::Digits2Raw() } else { Bool_t highGain = kFALSE; - // fill time bin values : - // 1. the signal starts at the time given by the digit - // 2. the rise is linear and the maximum is reached kTimePeak after start - // 3. the decay is gaussian with a sigma of kTimeRes - // 4. the signal is binned into kTimeBins bins - for (Int_t iTime = 0; iTime < kTimeBins; iTime++) { - Double_t time = iTime * kTimeMax/kTimeBins; - Int_t signal = 0; - if (time < digit->GetTime() + kTimePeak) { // signal is rising - signal = static_cast((fRan->Rndm() + digit->GetAmp()) * - (time - digit->GetTime() / kTimePeak) + 0.5); - } else { // signal is decaying - signal = static_cast((fRan->Rndm() + digit->GetAmp()) * - TMath::Gaus(time, digit->GetTime() + kTimePeak, kTimeRes) + 0.5); - } - if (signal < 0) - signal = 0; - adcValuesLow[iTime] = signal; - if (signal > 0x3FF) // larger than 10 bits - adcValuesLow[iTime] = 0x3FF; - adcValuesHigh[iTime] = static_cast(0.5 + (signal / kHighGainFactor)); - if (adcValuesHigh[iTime] > 0) - highGain = kTRUE; - } - // write low and eventually high gain channel - buffer->WriteChannel(idDDL, 0, 0, - kTimeBins, adcValuesLow, kThreshold); + buffer->WriteChannel(iDDL, 0, 0, + fkTimeBins, adcValuesLow, kThreshold); if (highGain) { - buffer->WriteChannel(idDDL, 0, kHighGainOffset, - kTimeBins, adcValuesHigh, 1); + buffer->WriteChannel(iDDL, 0, fHighGainOffset, + fkTimeBins, adcValuesHigh, 1); } } } - // write real header and close last file if (buffer) { buffer->Flush(); buffer->WriteDataHeader(kFALSE, kFALSE); delete buffer; } - gime->EmcalLoader()->UnloadDigits(); } @@ -310,6 +281,56 @@ AliLoader* AliEMCAL::MakeLoader(const char* topfoldername) return fLoader; } +//__________________________________________________________________ +Double_t AliEMCAL::RawResponseFunction(Double_t *x, Double_t *par) +{ + // Shape of the electronics raw reponse: + // 1. the signal rises linearly from par[4] to par[1] to reach the maximu par[3] + // 2. the signal decays with a gaussian shape for par[4]+par[1] with a sigma of par[2] + + Float_t xx = x[0] ; + + Double_t signal = 0. ; + + if (xx < par[4] + par[1]) // signal is rising + signal = (gRandom->Rndm() + par[3]) * (xx - par[4]) / (par[1] - par[4]) ; + else // signal is decaying + signal = (gRandom->Rndm() + par[3]) * TMath::Gaus(xx, par[4] + par[1], par[2]) ; + + return signal < 0. ? 0. : signal ; +} + +//__________________________________________________________________ +Bool_t AliEMCAL::RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const +{ + // for a start time dtime and an amplitude damp given by digit, + // calculates the raw sampled response + + const Int_t kRawSignalOverflow = 0x3FF ; + Bool_t highGain = kFALSE ; + + TF1 f1("signal", RawResponseFunction, 0, fkTimeBins, 5); + + f1.SetParNames("Time Max", "Peaking time", "Decay width", "Amp max", "Start time") ; + f1.SetParameter(0, fTimeMax) ; + f1.SetParameter(1, fTimePeak) ; + f1.SetParameter(2, fTimeRes) ; + f1.SetParameter(3, damp) ; + f1.SetParameter(4, dtime) ; + + for (Int_t iTime = 0; iTime < fkTimeBins; iTime++) { + Double_t time = iTime * fTimeMax/fkTimeBins; + Double_t signal = f1.Eval(time) ; + adcL[iTime] = static_cast(signal + 0.5) ; + if ( adcL[iTime] > kRawSignalOverflow) // larger than 10 bits + adcL[iTime] = kRawSignalOverflow ; + adcH[iTime] = static_cast(0.5 + (signal / fHighGainFactor)) ; + if (adcH[iTime] > 0) + highGain = kTRUE; + } + return highGain ; +} + //____________________________________________________________________________ void AliEMCAL::SetTreeAddress() { diff --git a/EMCAL/AliEMCAL.h b/EMCAL/AliEMCAL.h index ab642fc71c4..7e70cdd121e 100644 --- a/EMCAL/AliEMCAL.h +++ b/EMCAL/AliEMCAL.h @@ -45,21 +45,32 @@ class AliEMCAL : public AliDetector { {return AliEMCALGeometry::GetInstance(GetTitle(),"") ; } virtual void Hits2SDigits(); virtual Int_t IsVersion(void) const = 0 ; + Int_t GetRawFormatHighGainFactor() const { return fHighGainFactor ; } + Int_t GetRawFormatHighGainOffset() const { return fHighGainOffset ; } + Int_t GetRawFormatTimeBins() const { return fkTimeBins ; } + Double_t GetRawFormatTimeMax() const { return fTimeMax ; } + Double_t GetRawFormatTimePeak() const { return fTimePeak ; } + Double_t GetRawFormatTimeRes() const { return fTimeRes ; } virtual AliLoader* MakeLoader(const char* topfoldername); + static Double_t RawResponseFunction(Double_t *x, Double_t *par) ; virtual void SetTreeAddress() ; virtual const TString Version() const {return TString(" ") ; } AliEMCAL & operator = (const AliEMCAL & /*rvalue*/) { Fatal("operator =", "not implemented") ; return *this ; } - - protected: - //AliEMCALGeometry * fGeom ; // the geometry object + Bool_t RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const ; + Int_t fBirkC0; // constants for Birk's Law implementation Double_t fBirkC1; // constants for Birk's Law implementation Double_t fBirkC2; // constants for Birk's Law implementation - TRandom * fRan ; //! random number generator + Int_t fHighGainFactor ; // High gain attenuation factor of the raw RO signal + Int_t fHighGainOffset ; // offset added to the module id to distinguish high and low gain data + static const Int_t fkTimeBins = 256 ; // number of sampling bins of the raw RO signal + Double_t fTimeMax ; // maximum sampled time of the raw RO signal + Double_t fTimePeak ; // peaking time of the raw RO signal + Double_t fTimeRes ; // decay rime width of the raw RO signal ClassDef(AliEMCAL,7) // Electromagnetic calorimeter (base class) diff --git a/PHOS/AliPHOS.cxx b/PHOS/AliPHOS.cxx index b72a77ab934..42af88070d9 100644 --- a/PHOS/AliPHOS.cxx +++ b/PHOS/AliPHOS.cxx @@ -30,6 +30,8 @@ class TFile; #include #include #include +#include +#include #include // --- Standard library --- @@ -46,31 +48,33 @@ class TFile; ClassImp(AliPHOS) //____________________________________________________________________________ -AliPHOS:: AliPHOS() : AliDetector() + AliPHOS:: AliPHOS() : AliDetector() { // Default ctor - fName="PHOS"; + fName = "PHOS" ; fQATask = 0; fTreeQA = 0; fDebug = 0; - fRan = new TRandom(123456) ; } //____________________________________________________________________________ -AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title) +AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title) { // ctor : title is used to identify the layout - fQATask = 0; - fTreeQA = 0; - fDebug = 0; - fRan = new TRandom(123456) ; + fQATask = 0 ; + fTreeQA = 0 ; + fDebug = 0 ; + fTimeMax = 1.28E-5 ; + fTimePeak = 2.0E-6 ; + fTimeRes = 1.5E-6 ; + fHighGainFactor = 40; + fHighGainOffset = 0x200 ; } //____________________________________________________________________________ AliPHOS::~AliPHOS() { - delete fRan ; } //____________________________________________________________________________ @@ -409,19 +413,13 @@ void AliPHOS::Digits2Raw() } // some digitization constants - const Int_t kDDLOffset = 0x600; - const Double_t kTimeMax = 1.28E-5; - const Int_t kTimeBins = 256; - const Double_t kTimePeak = 2.0E-6; - const Double_t kTimeRes = 1.5E-6; - const Int_t kThreshold = 3; - const Int_t kHighGainFactor = 40; - const Int_t kHighGainOffset = 0x200; + const Int_t kDDLOffset = 0x600; // assigned to PHOS + const Int_t kThreshold = 3; // skip digits below this threshold AliAltroBuffer* buffer = NULL; Int_t prevDDL = -1; - Int_t adcValuesLow[kTimeBins]; - Int_t adcValuesHigh[kTimeBins]; + Int_t adcValuesLow[fkTimeBins]; + Int_t adcValuesHigh[fkTimeBins]; // loop over digits (assume ordered digits) for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { @@ -460,58 +458,34 @@ void AliPHOS::Digits2Raw() } // out of time range signal (?) - if (digit->GetTime() > kTimeMax) { + if (digit->GetTimeR() > fTimeMax) { buffer->FillBuffer(digit->GetAmp()); - buffer->FillBuffer(kTimeBins); // time bin + buffer->FillBuffer(fkTimeBins); // time bin buffer->FillBuffer(3); // bunch length buffer->WriteTrailer(3, relId[3], relId[2], module); // trailer // simulate linear rise and gaussian decay of signal } else { - Bool_t highGain = kFALSE; - - // fill time bin values : - // 1. the signal starts at the time given by the digit - // 2. the rise is linear and the maximum is reached kTimePeak after start - // 3. the decay is gaussian with a sigma of kTimeRes - // 4. the signal is binned into kTimeBins bins - for (Int_t iTime = 0; iTime < kTimeBins; iTime++) { - Double_t time = iTime * kTimeMax/kTimeBins; - Int_t signal = 0; - if (time < digit->GetTime() + kTimePeak) { // signal is rising - signal = static_cast(fRan->Rndm() + digit->GetAmp() * - (time - digit->GetTime()) / kTimePeak); - } else { // signal is decaying - signal = static_cast(fRan->Rndm() + digit->GetAmp() * - TMath::Gaus(time, digit->GetTime() + kTimePeak, kTimeRes)); - } - if (signal < 0) - signal = 0; - adcValuesLow[iTime] = signal; - if (signal > 0x3FF) // larger than 10 bits - adcValuesLow[iTime] = 0x3FF; - adcValuesHigh[iTime] = signal / kHighGainFactor; - if (adcValuesHigh[iTime] > 0) - highGain = kTRUE; - } - + Bool_t highGain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; + + // write low and eventually high gain channel buffer->WriteChannel(relId[3], relId[2], module, - kTimeBins, adcValuesLow, kThreshold); + fkTimeBins, adcValuesLow, kThreshold); if (highGain) { - buffer->WriteChannel(relId[3], relId[2], module + kHighGainOffset, - kTimeBins, adcValuesHigh, 1); + buffer->WriteChannel(relId[3], relId[2], module + fHighGainOffset, + fkTimeBins, adcValuesHigh, 1); } } } - + // write real header and close last file if (buffer) { buffer->Flush(); buffer->WriteDataHeader(kFALSE, kFALSE); delete buffer; } - + gime->PhosLoader()->UnloadDigits(); } @@ -535,6 +509,56 @@ AliLoader* AliPHOS::MakeLoader(const char* topfoldername) return fLoader; } +//__________________________________________________________________ +Double_t AliPHOS::RawResponseFunction(Double_t *x, Double_t *par) +{ + // Shape of the electronics raw reponse: + // 1. the signal rises linearly from par[4] to par[1] to reach the maximu par[3] + // 2. the signal decays with a gaussian shape for par[4]+par[1] with a sigma of par[2] + + Float_t xx = x[0] ; + + Double_t signal = 0. ; + + if (xx < par[4] + par[1]) // signal is rising + signal = (gRandom->Rndm() + par[3]) * (xx - par[4]) / (par[1] - par[4]) ; + else // signal is decaying + signal = (gRandom->Rndm() + par[3]) * TMath::Gaus(xx, par[4] + par[1], par[2]) ; + + return signal < 0. ? 0. : signal ; +} + +//__________________________________________________________________ +Bool_t AliPHOS::RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const +{ + // for a start time dtime and an amplitude damp given by digit, + // calculates the raw sampled response + + const Int_t kRawSignalOverflow = 0x3FF ; + Bool_t highGain = kFALSE ; + + TF1 f1("signal", RawResponseFunction, 0, fkTimeBins, 5); + + f1.SetParNames("Time Max", "Peaking time", "Decay width", "Amp max", "Start time") ; + f1.SetParameter(0, fTimeMax) ; + f1.SetParameter(1, fTimePeak) ; + f1.SetParameter(2, fTimeRes) ; + f1.SetParameter(3, damp) ; + f1.SetParameter(4, dtime) ; + + for (Int_t iTime = 0; iTime < fkTimeBins; iTime++) { + Double_t time = iTime * fTimeMax/fkTimeBins; + Double_t signal = f1.Eval(time) ; + adcL[iTime] = static_cast(signal + 0.5) ; + if ( adcL[iTime] > kRawSignalOverflow) // larger than 10 bits + adcL[iTime] = kRawSignalOverflow ; + adcH[iTime] = static_cast(0.5 + (signal / fHighGainFactor)) ; + if (adcH[iTime] > 0) + highGain = kTRUE; + } + return highGain ; +} + //____________________________________________________________________________ void AliPHOS::SetTreeAddress() { diff --git a/PHOS/AliPHOS.h b/PHOS/AliPHOS.h index a2b89d80ac5..0903baaaa87 100644 --- a/PHOS/AliPHOS.h +++ b/PHOS/AliPHOS.h @@ -48,8 +48,15 @@ public: {return AliPHOSGeometry::GetInstance(GetTitle(),"") ; } virtual void Hits2SDigits(); virtual Int_t IsVersion(void) const = 0 ; + Int_t GetRawFormatHighGainFactor() const { return fHighGainFactor ; } + Int_t GetRawFormatHighGainOffset() const { return fHighGainOffset ; } + Int_t GetRawFormatTimeBins() const { return fkTimeBins ; } + Double_t GetRawFormatTimeMax() const { return fTimeMax ; } + Double_t GetRawFormatTimePeak() const { return fTimePeak ; } + Double_t GetRawFormatTimeRes() const { return fTimeRes ; } virtual AliLoader* MakeLoader(const char* topfoldername); AliPHOSQAChecker * QAChecker() {return fQATask;} + static Double_t RawResponseFunction(Double_t *x, Double_t *par) ; virtual void SetTreeAddress(); virtual TTree * TreeQA() const {return fTreeQA; } virtual const TString Version() const {return TString(" ") ; } @@ -57,13 +64,22 @@ public: AliPHOS & operator = (const AliPHOS & /*rvalue*/) { Fatal("operator =", "not implemented") ; return *this ; } + protected: - + + Bool_t RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const ; + + AliPHOSQAChecker * fQATask ; //! PHOS checkers container TTree * fTreeQA ; // the QA tree that contains the alarms - TRandom * fRan ; //! random number generator + Int_t fHighGainFactor ; // High gain attenuation factor of the raw RO signal + Int_t fHighGainOffset ; // offset added to the module id to distinguish high and low gain data + static const Int_t fkTimeBins = 256 ; // number of sampling bins of the raw RO signal + Double_t fTimeMax ; // maximum sampled time of the raw RO signal + Double_t fTimePeak ; // peaking time of the raw RO signal + Double_t fTimeRes ; // decay rime width of the raw RO signal - ClassDef(AliPHOS,2) // Photon Spectrometer Detector (base class) + ClassDef(AliPHOS,3) // Photon Spectrometer Detector (base class) } ;