]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOS.cxx
Change for a realistic time distribution of the raw signal
[u/mrichter/AliRoot.git] / PHOS / AliPHOS.cxx
index fa2d48e40f586bde09abec7fc6cb1afea6473ab9..9eff25a38469336a04a9020ad364621ee40be6d9 100644 (file)
 
 // --- ROOT system ---
 class TFile;
-#include "TROOT.h"
-#include "TTree.h"
-#include "TFolder.h" 
+#include <TFolder.h> 
+#include <TTree.h>
+#include <TVirtualMC.h> 
+#include <TH1F.h> 
+#include <TF1.h> 
+#include <TRandom.h> 
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
+#include "AliMagF.h"
 #include "AliPHOS.h"
+#include "AliPHOSGetter.h"
 #include "AliRun.h"
-#include "AliMagF.h"
-#include "AliPHOSGeometry.h"
-#include "AliPHOSQAChecker.h" 
+#include "AliPHOSDigitizer.h"
+#include "AliPHOSSDigitizer.h"
+#include "AliPHOSDigit.h"
+#include "AliAltroBuffer.h"
+#include "AliLog.h"
 
 ClassImp(AliPHOS)
+
+Double_t AliPHOS::fgCapa        = 1.;        // 1pF 
+Int_t    AliPHOS::fgOrder       = 2 ;
+Double_t AliPHOS::fgTimeMax     = 2.56E-5 ;  // each sample is over 100 ns fTimeMax/fTimeBins
+Double_t AliPHOS::fgTimePeak    = 4.1E-6 ;   // 4 micro seconds
+Double_t AliPHOS::fgTimeTrigger = 100E-9 ;      // 100ns, just for a reference
+
 //____________________________________________________________________________
-AliPHOS:: AliPHOS() : AliDetector()
+  AliPHOS:: AliPHOS() : AliDetector()
 {
   // Default ctor
-  fName="PHOS";
+  fName   = "PHOS" ;
   fQATask = 0;
   fTreeQA = 0;
+  fDebug  = 0; 
+
 }
 
 //____________________________________________________________________________
-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;
+  fQATask = 0 ; 
+  fTreeQA = 0 ;
+  fDebug  = 0 ; 
+  fHighCharge        = 8.2 ;          // adjusted for a high gain range of 5.12 GeV (10 bits)
+  fHighGain          = 6.64 ; 
+  fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits) 
+  fLowGainOffset     = GetGeometry()->GetNModules() + 1 ;   
+    // offset added to the module id to distinguish high and low gain data
 }
 
 //____________________________________________________________________________
 AliPHOS::~AliPHOS() 
 {  
-  
+}
+
+//____________________________________________________________________________
+void AliPHOS::Copy(AliPHOS & phos)
+{
+  // copy method to be used byy the cpy ctor
+  TObject::Copy(phos) ; 
+  //  fQATask = AliPHOSQAChecker::Copy(*(phos.fQATask)) ; 
+  phos.fTreeQA = fTreeQA->CloneTree() ; 
+}
+
+//____________________________________________________________________________
+AliDigitizer* AliPHOS::CreateDigitizer(AliRunDigitizer* manager) const
+{
+  return new AliPHOSDigitizer(manager);
 }
 
 //____________________________________________________________________________
@@ -181,7 +218,7 @@ void AliPHOS::CreateMaterials()
   AliMaterial(15, "Ar$", 39.948, 18.0, dAr, 14.0, 0., 0, 0) ;   
  
  // ArCo2
-  Char_t namate[21];
+  Char_t namate[21]="";
   Float_t aGM[2] ; 
   Float_t zGM[2] ; 
   Float_t wGM[2] ; 
@@ -203,6 +240,7 @@ void AliPHOS::CreateMaterials()
   wGM[1] = 1. - arContent ;
   dGM    = wGM[0] * dAr + wGM[1] * dCO;
 
+  
   AliMixture(16, "ArCO2$", aGM, zGM, dGM,  2, wGM) ;
 
   // --- Stainless steel (let it be pure iron) ---
@@ -228,12 +266,14 @@ void AliPHOS::CreateMaterials()
   AliMixture(19, "Cables  $", aCA, zCA, dCA, -4, wCA) ;
 
 
-
   // --- Air ---
-  AliMaterial(99, "Air$", 14.61, 7.3, 0.001205, 30420., 67500., 0, 0) ;
-  
+  Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
+  Float_t zAir[4]={6.,7.,8.,18.};
+  Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
+  Float_t dAir = 1.20479E-3;
  
+  AliMixture(99, "Air$", aAir, zAir, dAir, 4, wAir) ;
+
   // DEFINITION OF THE TRACKING MEDIA
 
   // for PHOS: idtmed[699->798] equivalent to fIdtmed[0->100]
@@ -357,27 +397,216 @@ void AliPHOS::CreateMaterials()
 }
 
 //____________________________________________________________________________
-AliPHOSGeometry * AliPHOS::GetGeometry() const 
-{  
-  // gets the pointer to the AliPHOSGeometry unique instance 
+void AliPHOS::Digits2Raw()
+{
+// convert digits of the current event to raw data
+  
+  AliPHOSLoader * loader = dynamic_cast<AliPHOSLoader*>(fLoader) ; 
+
+  // get the digits
+  loader->LoadDigits();
+  TClonesArray* digits = loader->Digits() ;
+
+  if (!digits) {
+    Error("Digits2Raw", "no digits found !");
+    return;
+  }
+
+  // get the digitizer 
+  loader->LoadDigitizer();
+  AliPHOSDigitizer * digitizer = dynamic_cast<AliPHOSDigitizer *>(loader->Digitizer())  ; 
+  
+  // get the geometry
+  AliPHOSGeometry* geom = GetGeometry();
+  if (!geom) {
+    Error("Digits2Raw", "no geometry found !");
+    return;
+  }
+
+  // some digitization constants
+  const Int_t    kDDLOffset = 0x600; // assigned to PHOS
+  const Int_t    kThreshold = 1; // skip digits below this threshold
+
+  AliAltroBuffer* buffer = NULL;
+  Int_t prevDDL = -1;
+  Int_t adcValuesLow[fkTimeBins];
+  Int_t adcValuesHigh[fkTimeBins];
+
+  // loop over digits (assume ordered digits)
+  for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
+    AliPHOSDigit* digit = dynamic_cast<AliPHOSDigit *>(digits->At(iDigit)) ;
+    if (digit->GetAmp() < kThreshold) 
+      continue;
+    Int_t relId[4];
+    geom->AbsToRelNumbering(digit->GetId(), relId);
+    Int_t module = relId[0];
+   // Begin FIXME 
+    if (relId[1] != 0) 
+      continue;    // ignore digits from CPV
+   // End FIXME 
+
+    // PHOS EMCA has 4 DDL per module. Splitting is done based on the row number
+    Int_t iDDL = 4 * (module - 1) + (4 * (relId[2] - 1)) / geom->GetNPhi();
+
+    // new DDL
+    if (iDDL != prevDDL) {
+      // write real header and close previous file
+      if (buffer) {
+       buffer->Flush();
+       buffer->WriteDataHeader(kFALSE, kFALSE);
+       delete buffer;
+      }
+
+      // open new file and write dummy header
+      TString fileName("PHOS_") ;
+      fileName += (iDDL + kDDLOffset) ; 
+      fileName += ".ddl" ; 
+      buffer = new AliAltroBuffer(fileName.Data(), 1);
+      buffer->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
+
+      prevDDL = iDDL;
+    }
+
+    // out of time range signal (?)
+    if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
+      buffer->FillBuffer(digit->GetAmp());
+      buffer->FillBuffer(GetRawFormatTimeBins() );  // time bin
+      buffer->FillBuffer(3);          // bunch length      
+      buffer->WriteTrailer(3, relId[3], relId[2], module);  // trailer
+      
+    // calculate the time response function
+    } else {
+      Double_t energy = 0 ;  
+      if ( digit->GetId() <= geom->GetNModules() *  geom->GetNCristalsInModule())
+       energy = digit->GetAmp() * digitizer->GetEMCchannel() + digitizer->GetEMCpedestal() ; 
+      else 
+       energy = digit->GetAmp() * digitizer->GetCPVchannel() + digitizer->GetCPVpedestal() ;
+        
+      RawSampledResponse(digit->GetTimeR(), energy, adcValuesHigh, adcValuesLow) ; 
+      
+      buffer->WriteChannel(relId[3], relId[2], module, 
+                          GetRawFormatTimeBins(), adcValuesHigh, kThreshold);
+      buffer->WriteChannel(relId[3], relId[2], module + fLowGainOffset, 
+                          GetRawFormatTimeBins(), adcValuesLow, kThreshold);
+    }
+  }
+  
+  // write real header and close last file
+  if (buffer) {
+    buffer->Flush();
+    buffer->WriteDataHeader(kFALSE, kFALSE);
+    delete buffer;
+  }
   
-  return AliPHOSGeometry::GetInstance(GetTitle(),"") ;  
+  loader->UnloadDigits();
+}
+
+//____________________________________________________________________________
+void AliPHOS::Hits2SDigits()  
+{ 
+// create summable digits
 
+  AliPHOSSDigitizer* phosDigitizer = 
+    new AliPHOSSDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ;
+  phosDigitizer->SetEventRange(0, -1) ; // do all the events
+  phosDigitizer->ExecuteTask("all") ;
+}
+
+//____________________________________________________________________________
+AliLoader* AliPHOS::MakeLoader(const char* topfoldername)
+{
+//different behaviour than standard (singleton getter)
+// --> to be discussed and made eventually coherent
+ fLoader = new AliPHOSLoader(GetName(),topfoldername);
+ return fLoader;
+}
+
+//__________________________________________________________________
+Double_t AliPHOS::RawResponseFunction(Double_t *x, Double_t *par) 
+{
+  // Shape of the electronics raw reponse:
+  // It is a semi-gaussian, 2nd order Gamma function of the general form
+  // v(t) = n**n * Q * A**n / C *(t/tp)**n * exp(-n * t/tp) with 
+  // tp : peaking time par[0]
+  // n  : order of the function
+  // C  : integrating capacitor in the preamplifier
+  // A  : open loop gain of the preamplifier
+  // Q  : the total APD charge to be measured Q = C * energy
+  
+  Double_t signal ;
+  Double_t xx = x[0] - ( fgTimeTrigger + par[3] ) ; 
+
+  if (xx < 0 || xx > fgTimeMax) 
+    signal = 0. ;  
+  else { 
+    Double_t fac = par[0] * TMath::Power(fgOrder, fgOrder) * TMath::Power(par[1], fgOrder) / fgCapa ; 
+    signal = fac * par[2] * TMath::Power(xx / fgTimePeak, fgOrder) * TMath::Exp(-fgOrder * (xx / fgTimePeak)) ; 
+  }
+  return signal ;  
+}
+
+//__________________________________________________________________
+Double_t AliPHOS::RawResponseFunctionMax(Double_t charge, Double_t gain) 
+{
+  return ( charge * TMath::Power(fgOrder, fgOrder) * TMath::Power(gain, fgOrder) 
+     / ( fgCapa * TMath::Exp(fgOrder) ) );  
+
+}
+
+//__________________________________________________________________
+Bool_t AliPHOS::RawSampledResponse(const Double_t dtime, const Double_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 AliPHOS::RawResponseFunction
+
+  const Int_t kRawSignalOverflow = 0x3FF ; 
+  Bool_t lowGain = kFALSE ; 
+
+  TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
+
+  for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
+    signalF.SetParameter(0, GetRawFormatHighCharge() ) ; 
+    signalF.SetParameter(1, GetRawFormatHighGain() ) ; 
+    signalF.SetParameter(2, damp) ; 
+    signalF.SetParameter(3, dtime) ; 
+    Double_t time = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ;
+    Double_t signal = signalF.Eval(time) ;     
+    if ( static_cast<Int_t>(signal+0.5) > kRawSignalOverflow ){  // larger than 10 bits 
+      signal = kRawSignalOverflow ;
+      lowGain = kTRUE ; 
+    }
+    adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
+
+    signalF.SetParameter(0, GetRawFormatLowCharge() ) ;     
+    signalF.SetParameter(1, GetRawFormatLowGain() ) ; 
+    signal = signalF.Eval(time) ;  
+    if ( static_cast<Int_t>(signal+0.5) > kRawSignalOverflow)  // larger than 10 bits 
+      signal = kRawSignalOverflow ;
+    adcL[iTime] = static_cast<Int_t>(0.5 + signal ) ; 
+
+  }
+  return lowGain ; 
 }
 
 //____________________________________________________________________________
 void AliPHOS::SetTreeAddress()
 { 
-
+  // Links Hits in the Tree to Hits array
   TBranch *branch;
   char branchname[20];
   sprintf(branchname,"%s",GetName());
   
   // Branch address for hit tree
-  TTree *treeH = gAlice->TreeH();
-  if (treeH && fHits) {
+  TTree *treeH = TreeH();
+  if (treeH) {
     branch = treeH->GetBranch(branchname);
-    if (branch) branch->SetAddress(&fHits);
+    if (branch) 
+     { 
+       if (fHits == 0x0) fHits= new TClonesArray("AliPHOSHit",1000);
+       //Info("SetTreeAddress","<%s> Setting Hits Address",GetName());
+       branch->SetAddress(&fHits);
+     }
   }
 }
 
@@ -386,13 +615,24 @@ void AliPHOS::WriteQA()
 {
 
   // Make TreeQA in the output file. 
+
   if(fTreeQA == 0)
     fTreeQA = new TTree("TreeQA", "QA Alarms") ;    
   // Create Alarms branches
   Int_t bufferSize = 32000 ;    
   Int_t splitlevel = 0 ; 
-  TFolder * alarmsF = (TFolder*)gROOT->FindObjectAny("Folders/Run/Conditions/QA/PHOS") ; 
-  TString branchName(alarmsF->GetName());  
+
+  TFolder* topfold = GetLoader()->GetTopFolder(); //get top aliroot folder; skowron
+  TString phosqafn(AliConfig::Instance()->GetQAFolderName()+"/"); //get name of QAaut folder relative to top event; skowron
+  phosqafn+=GetName(); //hard wired string!!! add the detector name to the pathname; skowron 
+  TFolder * alarmsF = (TFolder*)topfold->FindObjectAny(phosqafn); //get the folder
+  if (alarmsF == 0x0)
+   {
+     Error("WriteQA","Can not find folder with qa alarms");
+     return;
+   }
+  TString branchName(alarmsF->GetName());
   TBranch * alarmsBranch = fTreeQA->Branch(branchName,"TFolder", &alarmsF, bufferSize, splitlevel);
   TString branchTitle = branchName + " QA alarms" ; 
   alarmsBranch->SetTitle(branchTitle);