]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOS.cxx
New AliTRDcluster constructor
[u/mrichter/AliRoot.git] / PHOS / AliPHOS.cxx
index bb76562e330ffa1d4c7fe549d65f3297f94ef51a..a5ee04d407c2c9103c9a0cede3749728593bfd51 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.111  2007/02/18 15:21:47  kharlov
+ * Corrections for memory leak in Digits2Raw due to AliAltroMapping
+ *
+ * Revision 1.110  2007/02/13 10:52:08  policheh
+ * Raw2SDigits() implemented
+ *
+ * Revision 1.109  2007/02/05 10:43:25  hristov
+ * Changes for correct initialization of Geant4 (Mihaela)
+ *
+ * Revision 1.108  2007/02/01 10:34:47  hristov
+ * Removing warnings on Solaris x86
+ *
+ * Revision 1.107  2007/01/29 16:29:37  kharlov
+ * Digits2Raw(): special workaround for digits with time out of range
+ *
+ * Revision 1.106  2007/01/17 17:28:56  kharlov
+ * Extract ALTRO sample generation to a separate class AliPHOSPulseGenerator
+ *
+ * Revision 1.105  2007/01/12 21:44:29  kharlov
+ * Simulate and reconstruct two gains simulaneouslsy
+ *
+ * Revision 1.104  2006/11/23 13:40:44  hristov
+ * Common class for raw data reading and ALTRO mappiing for PHOS and EMCAL (Gustavo, Cvetan)
+ *
+ * Revision 1.103  2006/11/14 17:11:15  hristov
+ * Removing inheritances from TAttLine, TAttMarker and AliRndm in AliModule. The copy constructor and assignment operators are moved to the private part of the class and not implemented. The corresponding changes are propagated to the detectors
+ *
+ * Revision 1.102  2006/10/27 17:14:27  kharlov
+ * Introduce AliDebug and AliLog (B.Polichtchouk)
+ *
+ * Revision 1.101  2006/10/13 06:47:29  kharlov
+ * Simulation of RAW data applies real mapping (B.Polichtchouk)
+ *
+ * Revision 1.100  2006/08/11 12:36:26  cvetan
+ * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
+ *
+ * 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 "AliPHOSDigitizer.h"
+#include "AliPHOSSDigitizer.h"
+#include "AliPHOSDigit.h"
+#include "AliAltroBuffer.h"
+#include "AliAltroMapping.h"
+#include "AliCaloAltroMapping.h"
+#include "AliLog.h"
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
+#include "AliCDBStorage.h"
+#include "AliPHOSCalibData.h"
+#include "AliPHOSPulseGenerator.h"
+#include "AliDAQ.h"
+#include "AliPHOSRawDecoder.h"
+#include "AliPHOSRawDigiProducer.h"
 
 ClassImp(AliPHOS)
 
 //____________________________________________________________________________
-AliPHOS::~AliPHOS()
+  AliPHOS:: AliPHOS() : AliDetector()
+{
+  // Default ctor
+  fName   = "PHOS" ;
+
+}
+
+//____________________________________________________________________________
+AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title)
+{
+  //   ctor : title is used to identify the layout
+}
+
+//____________________________________________________________________________
+AliPHOS::~AliPHOS() 
+{  
+}
+
+//____________________________________________________________________________
+AliDigitizer* AliPHOS::CreateDigitizer(AliRunDigitizer* manager) const
 {
-  // dtor
-  delete fHits ;
-  delete fDigits ;
+  return new AliPHOSDigitizer(manager);
 }
 
 //____________________________________________________________________________
@@ -105,11 +217,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} ;
@@ -118,12 +230,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.} ; 
@@ -160,40 +272,51 @@ 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]
-  Int_t * idtmed = fIdtmed->GetArray() - 699 ; 
   Int_t   isxfld = gAlice->Field()->Integ() ;
   Float_t sxmgmx = gAlice->Field()->Max() ;
 
@@ -261,11 +384,31 @@ 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) ;
+}
 
+//_____________________________________________________________________________
+void AliPHOS::Init()
+{
+  //
+  // Initialises cuts for PHOS
+  //
   // --- Set decent energy thresholds for gamma and electron tracking
+  Int_t * idtmed = fIdtmed->GetArray() - 699 ; 
 
   // Tracking threshold for photons and electrons in the scintillator crystal 
   gMC->Gstpar(idtmed[699], "CUTGAM",0.5E-4) ; 
@@ -297,31 +440,231 @@ void AliPHOS::CreateMaterials()
   gMC->Gstpar(idtmed[715], "LOSS",2.) ;
   gMC->Gstpar(idtmed[715], "DRAY",0.) ;
   gMC->Gstpar(idtmed[715], "STRA",2.) ;
+}
+
+//____________________________________________________________________________
+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 Float_t    kThreshold = 0.001; // skip digits below 1 MeV
+  const Int_t      kAdcThreshold = 1;  // Lower ADC threshold to write to raw data
+
+  Int_t prevDDL = -1;
+
+  // Create a shaper pulse object
+  AliPHOSPulseGenerator *pulse = new AliPHOSPulseGenerator();
+  
+  Int_t *adcValuesLow = new Int_t[pulse->GetRawFormatTimeBins()];
+  Int_t *adcValuesHigh= new Int_t[pulse->GetRawFormatTimeBins()];
+  
+  const Int_t maxDDL = 20;
+  AliAltroBuffer  *buffer[maxDDL];
+  AliAltroMapping *mapping[maxDDL];
+
+  for(Int_t jDDL=0; jDDL<maxDDL; jDDL++) {
+    buffer[jDDL]=0;
+    mapping[jDDL]=0;
+  }
+
+  //!!!!for debug!!!
+  Int_t modMax=-111;
+  Int_t colMax=-111;
+  Int_t rowMax=-111;
+  Float_t eMax=-333;
+  //!!!for debug!!!
+
+  // 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 
+
+    Int_t row = relId[2]-1;
+    Int_t col = relId[3]-1;
+
+    Int_t iRCU = -111;
+
+    //RCU0
+    if(0<=row&&row<32 && 0<=col&&col<28) iRCU=0;
+
+    //RCU1
+    if(0<=row&&row<32 && 28<=col&&col<56) iRCU=1;
+
+    //RCU2
+    if(32<=row&&row<64 && 0<=col&&col<28) iRCU=2;
+
+    //RCU3
+    if(32<=row&&row<64 && 28<=col&&col<56) iRCU=3;
+
+
+    // PHOS EMCA has 4 DDL per module. Splitting is based on the (row,column) numbers.
+    // PHOS internal convention: 1<module<5.
+    Int_t iDDL = 4 * (module - 1) + iRCU;
+
+    // new DDL
+    if (iDDL != prevDDL) {
+      if (buffer[iDDL] == 0) {
+       // open new file and write dummy header
+       TString fileName = AliDAQ::DdlFileName("PHOS",iDDL);
+
+       TString path = gSystem->Getenv("ALICE_ROOT");
+       path += "/PHOS/mapping/RCU";
+       path += iRCU;
+       path += ".data";
+
+       mapping[iDDL] = new AliCaloAltroMapping(path.Data());
+       buffer[iDDL]  = new AliAltroBuffer(fileName.Data(),mapping[iDDL]);
+       buffer[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
+      }
+      prevDDL = iDDL;
+    }
+
+    AliDebug(2,Form("digit E=%.4f GeV, t=%g s, (mod,col,row)=(%d,%d,%d)\n",
+                   digit->GetEnergy(),digit->GetTimeR(),
+                   relId[0]-1,relId[3]-1,relId[2]-1));
+    // if a signal is out of time range, write only trailer
+    if (digit->GetTimeR() > pulse->GetRawFormatTimeMax()*0.5 ) {
+      AliInfo("Signal is out of time range.\n");
+      buffer[iDDL]->FillBuffer(0);
+      buffer[iDDL]->FillBuffer(pulse->GetRawFormatTimeBins() );  // time bin
+      buffer[iDDL]->FillBuffer(3);                               // bunch length
+      buffer[iDDL]->WriteTrailer(3, relId[3]-1, relId[2]-1, 0);  // trailer
+      
+    // calculate the time response function
+    } else {
+      Double_t energy = 0 ;
+      module = relId[0];
+      if (digit->GetId() <= geom->GetNModules() * geom->GetNCristalsInModule()) {
+       energy=digit->GetEnergy();
+       if(energy>eMax) {eMax=energy; modMax=module; colMax=col; rowMax=row;}
+      }
+      else {
+       energy = 0; // CPV raw data format is now know yet
+      }
+      pulse->SetAmplitude(energy);
+      pulse->SetTZero(digit->GetTimeR());
+      pulse->MakeSamples();
+      pulse->GetSamples(adcValuesHigh, adcValuesLow) ; 
+      buffer[iDDL]->WriteChannel(relId[3]-1, relId[2]-1, 0, 
+                          pulse->GetRawFormatTimeBins(), adcValuesLow , kAdcThreshold);
+      buffer[iDDL]->WriteChannel(relId[3]-1, relId[2]-1, 1, 
+                          pulse->GetRawFormatTimeBins(), adcValuesHigh, kAdcThreshold);
+    }
+  }
+  
+  // write real header and close last file
+  for (Int_t iDDL=0; iDDL<maxDDL; iDDL++) {
+    if (buffer[iDDL]) {
+      buffer[iDDL]->Flush();
+      buffer[iDDL]->WriteDataHeader(kFALSE, kFALSE);
+      delete buffer[iDDL];
+      if (mapping[iDDL]) delete mapping[iDDL];
+    }
+  }
+  
+  AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n",
+        modMax,colMax,rowMax,eMax));
+
+  delete pulse;
+  loader->UnloadDigits();
+}
+
+//____________________________________________________________________________
+void AliPHOS::Hits2SDigits()  
+{ 
+// 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;
 }
 
 //____________________________________________________________________________
 void AliPHOS::SetTreeAddress()
 { 
-  // 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) ;
-  // }
+  // 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);
+     }
+  }
 }
 
+//____________________________________________________________________________
+Bool_t AliPHOS::Raw2SDigits(AliRawReader* rawReader)
+{
 
+  AliPHOSLoader * loader = dynamic_cast<AliPHOSLoader*>(fLoader) ;
+
+  TTree * tree = 0 ;
+  tree = loader->TreeS() ;
+  if ( !tree ) {
+    loader->MakeTree("S");
+    tree = loader->TreeS() ;
+  }
+
+  TClonesArray * sdigits = loader->SDigits() ;
+  if(!sdigits) { 
+    loader->MakeSDigitsArray();
+    sdigits = loader->SDigits();
+  }
+  sdigits->Clear();
+
+  AliPHOSRawDecoder dc(rawReader);
+  AliPHOSRawDigiProducer pr;
+  pr.MakeDigits(sdigits,&dc);
+
+  Int_t bufferSize = 32000 ;
+  TBranch * sdigitsBranch = tree->Branch("PHOS",&sdigits,bufferSize);
+  tree->Fill();
+
+  fLoader->WriteSDigits("OVERWRITE");
+  return kTRUE;
+    
+}