]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOS.cxx
fixing wrong define
[u/mrichter/AliRoot.git] / PHOS / AliPHOS.cxx
index 45c7442b7cdfb31f4097fbb138963929c0a93793..8aecbb2211ce001cd7e2982a2a8ae63935ff03ac 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-
 /* $Id$ */
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.99  2006/06/28 11:36:09  cvetan
+ * New detector numbering scheme (common for DAQ/HLT/Offline). All the subdetectors shall use the AliDAQ class for the sim and rec of the raw data. The AliDAQ and raw reader classes now provide all the necessary interfaces to write and select the detector specific raw-data payload. Look into the AliDAQ.h and AliRawReader.h for more details.
+ *
+ * Revision 1.98  2006/05/11 11:30:48  cvetan
+ * Major changes in AliAltroBuffer. Now it can be used only for writing of raw data. All the corresponding read method are removed. It is based now on AliFstream in order to avoid endianess problems. The altro raw data is written always with little endian
+ *
+ * Revision 1.97  2006/04/22 10:30:17  hristov
+ * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
+ *
+ * Revision 1.96  2006/04/07 08:41:59  hristov
+ * Follow AliAlignObj framework and remove AliPHOSAlignData (Yu.Kharlov)
+ *
+ * Revision 1.95  2006/03/14 19:40:41  kharlov
+ * Remove De-digitizing of raw data and digitizing the raw data fit
+ *
+ * Revision 1.94  2006/03/07 18:56:25  kharlov
+ * CDB is passed via environment variable
+ *
+ * Revision 1.93  2005/11/22 08:45:11  kharlov
+ * Calibration is read from CDB if any (Boris Polichtchouk)
+ *
+ * Revision 1.92  2005/11/03 13:09:19  hristov
+ * Removing meaningless const declarations (linuxicc)
+ *
+ * Revision 1.91  2005/07/27 15:08:53  kharlov
+ * Mixture ArCO2 is corrected
+ *
+ * Revision 1.90  2005/06/17 07:39:07  hristov
+ * Removing GetDebug and SetDebug from AliRun and AliModule. Using AliLog for the messages
+ *
+ * Revision 1.89  2005/05/28 12:10:07  schutz
+ * Copy constructor is corrected (by T.P.)
+ *
+ */
 
 //_________________________________________________________________________
 // Base Class for PHOS description:
 //*-- Author: Laurent Aphecetche & Yves Schutz (SUBATECH) 
 //////////////////////////////////////////////////////////////////////////////
 
+
 // --- ROOT system ---
+class TFile;
+#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 "AliMC.h"
+#include "AliPHOSGetter.h"
 #include "AliRun.h"
-#include "AliMagF.h"
+#include "AliPHOSDigitizer.h"
+#include "AliPHOSSDigitizer.h"
+#include "AliPHOSDigit.h"
+#include "AliAltroBuffer.h"
+#include "AliLog.h"
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
+#include "AliCDBStorage.h"
+#include "AliPHOSCalibData.h"
+#include "AliDAQ.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
+
+Double_t AliPHOS::fgHighCharge  = 8.2;       // adjusted for a high gain range of 5.12 GeV (10 bits)
+Double_t AliPHOS::fgHighGain    = 6.64;
+Double_t AliPHOS::fgHighLowGainFactor = 16.; // adjusted for a low gain range of 82 GeV (10 bits) 
+
 //____________________________________________________________________________
-AliPHOS::AliPHOS():AliDetector()
+  AliPHOS:: AliPHOS() : AliDetector()
 {
-  // ctor 
-  //We do not create objects, because these pointers will be overwritten durin reading from file.
-  fEmcRecPoints  = 0 ; 
-  fPpsdRecPoints = 0 ;
-  fTrackSegments = 0 ;
-  fRecParticles  = 0 ;
+  // Default ctor
+  fName   = "PHOS" ;
 
 }
+
 //____________________________________________________________________________
-AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name,title) 
+AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title)
 {
-  // ctor
-  
-  fEmcRecPoints  = new AliPHOSRecPoint::RecPointsList(10) ; 
-  fPpsdRecPoints = new AliPHOSRecPoint::RecPointsList(10) ;
-  fTrackSegments = new AliPHOSTrackSegment::TrackSegmentsList("AliPHOSTrackSegment", 10) ;
-  fRecParticles  = new AliPHOSRecParticle::RecParticlesList("AliPHOSRecParticle", 10) ;
-  fDebugLevel    = 0;
-  
+  //   ctor : title is used to identify the layout
+}
+
+//____________________________________________________________________________
+AliPHOS::~AliPHOS() 
+{  
+}
+
+//____________________________________________________________________________
+void AliPHOS::Copy(TObject &obj)const
+{
+  // copy method to be used by the cpy ctor
+  TObject::Copy(obj);
 }
+
 //____________________________________________________________________________
-AliPHOS::~AliPHOS()
+AliDigitizer* AliPHOS::CreateDigitizer(AliRunDigitizer* manager) const
 {
-  // dtor
-  if(fEmcRecPoints)
-    fEmcRecPoints->Delete() ;
-  delete fEmcRecPoints ;
-  if(fPpsdRecPoints)
-    fPpsdRecPoints->Delete() ;
-  delete fPpsdRecPoints ;
-  if(fTrackSegments)
-    fTrackSegments->Delete() ;
-  delete fTrackSegments ;
-  if(fRecParticles)
-    fRecParticles->Delete() ;
-  delete fRecParticles ;
-  delete fHits ;
-  delete fDigits ;
+  return new AliPHOSDigitizer(manager);
 }
 
 //____________________________________________________________________________
@@ -141,11 +193,11 @@ void AliPHOS::CreateMaterials()
   Float_t aTI[2] = {12.011, 1.00794} ;
   Float_t zTI[2] = {6.0, 1.0} ;
   Float_t wTI[2] = {1.0, 1.0} ;
-  Float_t dTI = 0.1 ;
+  Float_t dTI = 0.04 ;
 
   AliMixture(7, "Thermo Insul.$", aTI, zTI, dTI, -2, wTI) ;
 
-  // --- Textolitn ---
+  // --- Textolith ---
   Float_t aTX[4] = {16.0, 28.09, 12.011, 1.00794} ;
   Float_t zTX[4] = {8.0, 14.0, 6.0, 1.0} ;
   Float_t wTX[4] = {292.0, 68.0, 462.0, 736.0} ;
@@ -154,12 +206,12 @@ void AliPHOS::CreateMaterials()
   AliMixture(8, "Textolit$", aTX, zTX, dTX, -4, wTX) ;
 
   //--- FR4  ---
-  Float_t aFR[3] = {28.0855, 15.9994, 17.749} ; 
-  Float_t zFR[3] = {14., 8., 8.875} ; 
-  Float_t wFR[3] = {.28, .32, .4} ;
+  Float_t aFR[4] = {16.0, 28.09, 12.011, 1.00794} ;
+  Float_t zFR[4] = {8.0, 14.0, 6.0, 1.0} ;
+  Float_t wFR[4] = {292.0, 68.0, 462.0, 736.0} ;
   Float_t dFR = 1.8 ; 
 
-  AliMixture(9, "FR4$", aFR, zFR, dFR, -3, wFR) ;
+  AliMixture(9, "FR4$", aFR, zFR, dFR, -4, wFR) ;
 
   // --- The Composite Material for  micromegas (so far polyetylene) ---                                       
   Float_t aCM[2] = {12.01, 1.} ; 
@@ -196,36 +248,48 @@ void AliPHOS::CreateMaterials()
   Float_t dAr = 0.001782 ; 
   AliMaterial(15, "Ar$", 39.948, 18.0, dAr, 14.0, 0., 0, 0) ;   
  
- // ArCo2
-  Char_t namate[21];
-  Float_t aGM[2] ; 
-  Float_t zGM[2] ; 
-  Float_t wGM[2] ; 
-  Float_t dGM ; 
+  // Ar+CO2 Mixture (80% / 20%)
+  Float_t arContent = 0.80 ;  // Ar-content of the ArCO2-mixture
+  Float_t aArCO[3]  = {39.948, 12.0, 16.0} ;
+  Float_t zArCO[3]  = {18.0  ,  6.0,  8.0} ;
+  Float_t wArCO[3];
+  wArCO[0] = arContent;
+  wArCO[1] = (1-arContent)*1;
+  wArCO[2] = (1-arContent)*2;
+  Float_t dArCO = arContent*dAr + (1-arContent)*dCO ;
+  AliMixture(16, "ArCO2$", aArCO, zArCO, dArCO,  -3, wArCO) ;
 
-  Float_t absL, radL, density ;
-  Float_t buf[1] ;
-  Int_t nbuf ;
+  // --- Stainless steel (let it be pure iron) ---
+  AliMaterial(17, "Steel$", 55.845, 26, 7.87, 1.76, 0., 0, 0) ;
 
-  gMC->Gfmate((*fIdmate)[15], namate, aGM[0], zGM[0], density, radL, absL, buf, nbuf) ; // Get properties of Ar 
-  gMC->Gfmate((*fIdmate)[14], namate, aGM[1], zGM[1], density, radL, absL, buf, nbuf) ; // Get properties of CO2 
 
+  // --- Fiberglass ---
+  Float_t aFG[4] = {16.0, 28.09, 12.011, 1.00794} ;
+  Float_t zFG[4] = {8.0, 14.0, 6.0, 1.0} ;
+  Float_t wFG[4] = {292.0, 68.0, 462.0, 736.0} ;
+  Float_t dFG    = 1.9 ;
 
-  // Create gas mixture 
+  AliMixture(18, "Fibergla$", aFG, zFG, dFG, -4, wFG) ;
 
-  Float_t arContent    = 0.80 ;  // Ar-content of the Ar/CO2-mixture (80% / 20%) 
-  wGM[0] = arContent;
-  wGM[1] = 1. - arContent ;
-  dGM    = wGM[0] * dAr + wGM[1] * dCO;
+  // --- Cables in Air box  ---
+  // SERVICES
+
+  Float_t aCA[4] = { 1.,12.,55.8,63.5 };
+  Float_t zCA[4] = { 1.,6.,26.,29. }; 
+  Float_t wCA[4] = { .014,.086,.42,.48 };
+  Float_t dCA    = 0.8 ;  //this density is raw estimation, if you know better - correct
+
+  AliMixture(19, "Cables  $", aCA, zCA, dCA, -4, wCA) ;
 
-  AliMixture(16, "ArCO2$", aGM, zGM, dGM,  2, wGM) ;
 
   // --- 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]
@@ -297,6 +361,18 @@ void AliPHOS::CreateMaterials()
   AliMedium(16, "ArCo2      $", 16, 1,
             isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.01, 0, 0) ;
  
+  // Stainless steel                                                                -> idtmed[716]
+  AliMedium(17, "Steel     $", 17, 0,
+            isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.0001, 0, 0) ;
+
+  // Fibergalss                                                                     -> idtmed[717]
+  AliMedium(18, "Fiberglass$", 18, 0,
+            isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ;
+
+  // Cables in air                                                                  -> idtmed[718]
+  AliMedium(19, "Cables    $", 19, 0,
+            isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ;
+
   // Air                                                                            -> idtmed[798] 
   AliMedium(99, "Air          $", 99, 0,
             isxfld, sxmgmx, 10.0, 1.0, 0.1, 0.1, 10.0, 0, 0) ;
@@ -337,27 +413,219 @@ void AliPHOS::CreateMaterials()
 }
 
 //____________________________________________________________________________
-void AliPHOS::SetTreeAddress()
+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) {
+    AliError(Form("No digits found !"));
+    return;
+  }
+
+  // get the geometry
+  AliPHOSGeometry* geom = GetGeometry();
+  if (!geom) {
+    AliError(Form("No geometry found !"));
+    return;
+  }
+
+  // some digitization constants
+//   const Int_t    kThreshold = 1; // skip digits below this threshold // YVK
+  const Float_t    kThreshold = 0.001; // skip digits below 1 MeV
+  const Int_t      kAdcThreshold = 1;  // Lower ADC threshold to write to raw data
+
+  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->GetEnergy() < 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 = AliDAQ::DdlFileName("PHOS",iDDL);
+      buffer = new AliAltroBuffer(fileName.Data());
+      buffer->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
+
+      prevDDL = iDDL;
+    }
+
+    // out of time range signal (?)
+    if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
+      buffer->FillBuffer((Int_t)digit->GetEnergy());
+      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 ;
+      Int_t   module = relId[0];
+      if ( digit->GetId() <= geom->GetNModules() *  geom->GetNCristalsInModule()) {
+       energy=digit->GetEnergy();
+      }
+      else {
+//     energy = digit->GetAmp()*digitizer->GetCPVchannel()+digitizer->GetCPVpedestal();
+       energy = 0; // CPV raw data format is now know yet
+      }        
+      Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), energy, adcValuesHigh, adcValuesLow) ; 
+      
+      if (lowgain) 
+       buffer->WriteChannel(relId[3], relId[2], module + GetGeometry()->GetNModules() + 1, 
+                            GetRawFormatTimeBins(), adcValuesLow , kAdcThreshold);
+      else 
+       buffer->WriteChannel(relId[3], relId[2], module, 
+                            GetRawFormatTimeBins(), adcValuesHigh, kAdcThreshold);
+      
+    }
+  }
+  
+  // write real header and close last file
+  if (buffer) {
+    buffer->Flush();
+    buffer->WriteDataHeader(kFALSE, kFALSE);
+    delete buffer;
+  }
+  
+  loader->UnloadDigits();
+}
+
+//____________________________________________________________________________
+void AliPHOS::Hits2SDigits()  
 { 
-  // TBranch *branch;
-  AliDetector::SetTreeAddress();
-
- //  //Branch address for TreeR: RecParticles
-//   TTree *treeR = gAlice->TreeR();
-//   if ( treeR && fRecParticles ) {
-//     branch = treeR->GetBranch("PHOSRP");
-//     if (branch) branch->SetAddress(&fRecParticles) ;
-//   }
-//    //Branch address for TreeR: TrackSegments
-//   if ( treeR && fTrackSegments ) {
-//     branch = treeR->GetBranch("PHOSTS");
-//     if (branch) branch->SetAddress(&fTrackSegments) ;
-//   }
-//   //Branch address for TreeR: EmcRecPoint
-//  if ( treeR && fEmcRecPoints ) {
-//     branch = treeR->GetBranch("PHOSEmcRP");
-//     if (branch) branch->SetAddress(&fEmcRecPoints) ;
-  // }
+// create summable digits
+
+  AliPHOSSDigitizer phosDigitizer(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(Double_t dtime, 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
+  // Input: dtime - signal start time
+  //        damp  - signal amplitude (energy)
+  // Output: adcH - array[fkTimeBins] of 10-bit samples for high-gain channel
+  //         adcL - array[fkTimeBins] of 10-bit samples for low-gain channel
+
+  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 = TreeH();
+  if (treeH) {
+    branch = treeH->GetBranch(branchname);
+    if (branch) 
+     { 
+       if (fHits == 0x0) fHits= new TClonesArray("AliPHOSHit",1000);
+       //AliInfo(Form("<%s> Setting Hits Address",GetName()));
+       branch->SetAddress(&fHits);
+     }
+  }
+}