]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOS.cxx
Implementation of copy constructor and assignement operators (Marian)
[u/mrichter/AliRoot.git] / PHOS / AliPHOS.cxx
index 2ea0699fe562443d8bbead5fda1226a755d6240f..d4131c882e7edf15082de71ffada55656da8e1c2 100644 (file)
 /* History of cvs commits:
  *
  * $Log$
- * 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.116  2007/10/10 09:05:10  schutz
+ * Changing name QualAss to QA
  *
- * 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.115  2007/08/22 09:20:50  hristov
+ * Updated QA classes (Yves)
  *
- * Revision 1.102  2006/10/27 17:14:27  kharlov
- * Introduce AliDebug and AliLog (B.Polichtchouk)
+ * Revision 1.114  2007/08/07 14:12:03  kharlov
+ * Quality assurance added (Yves Schutz)
  *
- * Revision 1.101  2006/10/13 06:47:29  kharlov
- * Simulation of RAW data applies real mapping (B.Polichtchouk)
+ * Revision 1.113  2007/07/18 16:29:54  policheh
+ * Raw Sdigits energy converted to GeV.
  *
- * 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.112  2007/02/25 22:59:13  policheh
+ * Digits2Raw(): ALTRO buffer and mapping created per each DDL.
  *
- * 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.111  2007/02/18 15:21:47  kharlov
+ * Corrections for memory leak in Digits2Raw due to AliAltroMapping
  *
- * 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.110  2007/02/13 10:52:08  policheh
+ * Raw2SDigits() implemented
  *
- * 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.109  2007/02/05 10:43:25  hristov
+ * Changes for correct initialization of Geant4 (Mihaela)
  *
- * Revision 1.96  2006/04/07 08:41:59  hristov
- * Follow AliAlignObj framework and remove AliPHOSAlignData (Yu.Kharlov)
+ * Revision 1.108  2007/02/01 10:34:47  hristov
+ * Removing warnings on Solaris x86
  *
- * Revision 1.95  2006/03/14 19:40:41  kharlov
- * Remove De-digitizing of raw data and digitizing the raw data fit
+ * Revision 1.107  2007/01/29 16:29:37  kharlov
+ * Digits2Raw(): special workaround for digits with time out of range
  *
- * 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.)
+ * 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
  */
 
 //_________________________________________________________________________
@@ -104,20 +91,16 @@ class TFile;
 #include "AliCDBEntry.h"
 #include "AliCDBStorage.h"
 #include "AliPHOSCalibData.h"
+#include "AliPHOSPulseGenerator.h"
 #include "AliDAQ.h"
+#include "AliPHOSRawDecoder.h"
+#include "AliPHOSRawDigiProducer.h"
+#include "AliPHOSQAChecker.h"
+#include "AliPHOSRecoParamEmc.h"
+#include "AliPHOSSimParam.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()
 {
@@ -238,7 +221,8 @@ void AliPHOS::CreateMaterials()
   Float_t aG10[4] = { 12., 1., 16., 28.} ;
   Float_t zG10[4] = { 6., 1., 8., 14.} ;
   Float_t wG10[4] = { .259, .288, .248, .205} ;
-  Float_t dG10  = 1.7 ;
+  Float_t dG10  = 1.7 ; 
+
   
   AliMixture(12, "G10$", aG10, zG10, dG10, -4, wG10);
 
@@ -267,7 +251,8 @@ void AliPHOS::CreateMaterials()
   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) ;
+  AliMixture(16, "ArCO2$", aArCO, zArCO, dArCO,  -3, wArCO) ; 
+
 
   // --- Stainless steel (let it be pure iron) ---
   AliMaterial(17, "Steel$", 55.845, 26, 7.87, 1.76, 0., 0, 0) ;
@@ -303,7 +288,6 @@ void AliPHOS::CreateMaterials()
   // 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() ;
 
@@ -386,8 +370,16 @@ void AliPHOS::CreateMaterials()
   // 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) ; 
@@ -419,7 +411,6 @@ void AliPHOS::CreateMaterials()
   gMC->Gstpar(idtmed[715], "LOSS",2.) ;
   gMC->Gstpar(idtmed[715], "DRAY",0.) ;
   gMC->Gstpar(idtmed[715], "STRA",2.) ;
-
 }
 
 //____________________________________________________________________________
@@ -427,6 +418,10 @@ void AliPHOS::Digits2Raw()
 {
 // convert digits of the current event to raw data
   
+  if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
+    AliError("Energy digitization should be OFF if use Digits2Raw") ;
+  }
+
   AliPHOSLoader * loader = dynamic_cast<AliPHOSLoader*>(fLoader) ; 
 
   // get the digits
@@ -445,16 +440,30 @@ void AliPHOS::Digits2Raw()
     return;
   }
 
+  // get mapping from OCDB
+  const TObjArray* maps = AliPHOSRecoParamEmc::GetMappings();
+  if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
+
   // 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];
 
+  // 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;
@@ -472,10 +481,10 @@ void AliPHOS::Digits2Raw()
     geom->AbsToRelNumbering(digit->GetId(), relId);
     Int_t module = relId[0];
  
-   // Begin FIXME 
+    // Begin FIXME 
     if (relId[1] != 0) 
       continue;    // ignore digits from CPV
-   // End FIXME 
+    // End FIXME 
 
     Int_t row = relId[2]-1;
     Int_t col = relId[3]-1;
@@ -501,68 +510,69 @@ void AliPHOS::Digits2Raw()
 
     // new DDL
     if (iDDL != prevDDL) {
-      // write real header and close previous file
-      if (buffer) {
-       buffer->Flush();
-       buffer->WriteDataHeader(kFALSE, kFALSE);
-       delete buffer;
+      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] = (AliAltroMapping*)maps->At(iRCU);
+       buffer[iDDL]  = new AliAltroBuffer(fileName.Data(),mapping[iDDL]);
+       buffer[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
       }
-
-      // 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";
-
-      AliAltroMapping* mapping = new AliCaloAltroMapping(path.Data());
-      buffer = new AliAltroBuffer(fileName.Data(),mapping);
-      buffer->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
-
       prevDDL = iDDL;
     }
 
-    // out of time range signal (?)
-    if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
+    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->FillBuffer((Int_t)digit->GetEnergy());
-      buffer->FillBuffer(GetRawFormatTimeBins() );  // time bin
-      buffer->FillBuffer(3);          // bunch length      
-      buffer->WriteTrailer(3, relId[3], relId[2], module);  // trailer
+      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 ;
-      Int_t   module = relId[0];
-      if ( digit->GetId() <= geom->GetNModules() *  geom->GetNCristalsInModule()) {
+      module = relId[0];
+      if (digit->GetId() <= geom->GetNModules() * geom->GetNCristalsInModule()) {
        energy=digit->GetEnergy();
-       AliDebug(2,Form("digit energy: %f\n",digit->GetEnergy()));
        if(energy>eMax) {eMax=energy; modMax=module; colMax=col; rowMax=row;}
       }
       else {
        energy = 0; // CPV raw data format is now know yet
-      }        
-      Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), energy, adcValuesHigh, adcValuesLow) ; 
-      
-       buffer->WriteChannel(relId[3]-1, relId[2]-1, 0, 
-                            GetRawFormatTimeBins(), adcValuesLow , kAdcThreshold);
-       buffer->WriteChannel(relId[3]-1, relId[2]-1, 1, 
-                            GetRawFormatTimeBins(), adcValuesHigh, kAdcThreshold);
-      
+      }
+      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
-  if (buffer) {
-    buffer->Flush();
-    buffer->WriteDataHeader(kFALSE, kFALSE);
-    delete buffer;
+  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();
 }
 
@@ -573,9 +583,11 @@ void AliPHOS::Hits2SDigits()
 
   AliPHOSSDigitizer phosDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ;
   phosDigitizer.SetEventRange(0, -1) ; // do all the events
   phosDigitizer.ExecuteTask("all") ; 
 }
 
+
 //____________________________________________________________________________
 AliLoader* AliPHOS::MakeLoader(const char* topfoldername)
 {
@@ -585,79 +597,6 @@ AliLoader* AliPHOS::MakeLoader(const char* 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) ;
-    AliDebug(4,Form("iTime: %d Energy: %f HG signal: %f adcH: %d ",iTime,damp,signal,adcH[iTime]));
-
-    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 ) ; 
-    AliDebug(4,Form("..LG: %f adcL: %d\n",signal,adcL[iTime]));
-
-  }
-  return lowGain ; 
-}
-
 //____________________________________________________________________________
 void AliPHOS::SetTreeAddress()
 { 
@@ -678,3 +617,36 @@ void AliPHOS::SetTreeAddress()
   }
 }
 
+//____________________________________________________________________________
+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->Branch("PHOS",&sdigits,bufferSize);
+  tree->Fill();
+
+  fLoader->WriteSDigits("OVERWRITE");
+  return kTRUE;
+    
+}