Made the raw data format conversion from digits more modular (+2 more methods) to...
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 16 Jul 2004 17:00:04 +0000 (17:00 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 16 Jul 2004 17:00:04 +0000 (17:00 +0000)
EMCAL/AliEMCAL.cxx
EMCAL/AliEMCAL.h
PHOS/AliPHOS.cxx
PHOS/AliPHOS.h

index 3b6dd1d..1e95f60 100644 (file)
@@ -30,6 +30,9 @@ class TFile;
 #include <TFolder.h> 
 #include <TTree.h>
 #include <TVirtualMC.h> 
+#include <TH1F.h> 
+#include <TF1.h> 
+#include <TRandom.h> 
 
 // --- 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<Int_t>((fRan->Rndm() + digit->GetAmp()) * 
-                                     (time - digit->GetTime() / kTimePeak) + 0.5);
-       } else {                                        // signal is decaying
-         signal = static_cast<Int_t>((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<Int_t>(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<Int_t>(signal + 0.5) ;
+    if ( adcL[iTime] > kRawSignalOverflow)  // larger than 10 bits 
+      adcL[iTime] = kRawSignalOverflow ;
+    adcH[iTime] = static_cast<Int_t>(0.5 + (signal / fHighGainFactor)) ;
+    if (adcH[iTime] > 0) 
+      highGain = kTRUE;  
+  }
+  return highGain ; 
+}
+
 //____________________________________________________________________________
 void AliEMCAL::SetTreeAddress()
 { 
index ab642fc..7e70cdd 100644 (file)
@@ -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)
 
index b72a77a..42af880 100644 (file)
@@ -30,6 +30,8 @@ class TFile;
 #include <TFolder.h> 
 #include <TTree.h>
 #include <TVirtualMC.h> 
+#include <TH1F.h> 
+#include <TF1.h> 
 #include <TRandom.h> 
 
 // --- 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<Int_t>(fRan->Rndm() + digit->GetAmp() * 
-                        (time - digit->GetTime()) / kTimePeak);
-       } else {                                        // signal is decaying
-         signal = static_cast<Int_t>(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<Int_t>(signal + 0.5) ;
+    if ( adcL[iTime] > kRawSignalOverflow)  // larger than 10 bits 
+      adcL[iTime] = kRawSignalOverflow ;
+    adcH[iTime] = static_cast<Int_t>(0.5 + (signal / fHighGainFactor)) ;
+    if (adcH[iTime] > 0) 
+      highGain = kTRUE;  
+  }
+  return highGain ; 
+}
+
 //____________________________________________________________________________
 void AliPHOS::SetTreeAddress()
 { 
index a2b89d8..0903baa 100644 (file)
@@ -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)
 
 } ;