]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOS.cxx
Coverity fix
[u/mrichter/AliRoot.git] / PHOS / AliPHOS.cxx
index 6556010ee37d957eb5280e440cbec3abf02d65b3..8c0e9e8fddc460ccdd7e0dbe4b28c7607744d97a 100644 (file)
 /* History of cvs commits:
  *
  * $Log$
- * 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.116  2007/10/10 09:05:10  schutz
+ * Changing name QualAss to QA
  *
- * 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.115  2007/08/22 09:20:50  hristov
+ * Updated QA classes (Yves)
  *
- * 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.114  2007/08/07 14:12:03  kharlov
+ * Quality assurance added (Yves Schutz)
  *
- * Revision 1.96  2006/04/07 08:41:59  hristov
- * Follow AliAlignObj framework and remove AliPHOSAlignData (Yu.Kharlov)
+ * Revision 1.113  2007/07/18 16:29:54  policheh
+ * Raw Sdigits energy converted to GeV.
  *
- * Revision 1.95  2006/03/14 19:40:41  kharlov
- * Remove De-digitizing of raw data and digitizing the raw data fit
+ * Revision 1.112  2007/02/25 22:59:13  policheh
+ * Digits2Raw(): ALTRO buffer and mapping created per each DDL.
  *
- * Revision 1.94  2006/03/07 18:56:25  kharlov
- * CDB is passed via environment variable
+ * Revision 1.111  2007/02/18 15:21:47  kharlov
+ * Corrections for memory leak in Digits2Raw due to AliAltroMapping
  *
- * Revision 1.93  2005/11/22 08:45:11  kharlov
- * Calibration is read from CDB if any (Boris Polichtchouk)
+ * Revision 1.110  2007/02/13 10:52:08  policheh
+ * Raw2SDigits() implemented
  *
- * Revision 1.92  2005/11/03 13:09:19  hristov
- * Removing meaningless const declarations (linuxicc)
+ * Revision 1.109  2007/02/05 10:43:25  hristov
+ * Changes for correct initialization of Geant4 (Mihaela)
  *
- * Revision 1.91  2005/07/27 15:08:53  kharlov
- * Mixture ArCO2 is corrected
+ * Revision 1.108  2007/02/01 10:34:47  hristov
+ * Removing warnings on Solaris x86
  *
- * 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.107  2007/01/29 16:29:37  kharlov
+ * Digits2Raw(): special workaround for digits with time out of range
  *
- * 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
  */
 
 //_________________________________________________________________________
 
 // --- ROOT system ---
 class TFile;
+#include <TF1.h> 
 #include <TFolder.h> 
-#include <TTree.h>
-#include <TVirtualMC.h> 
+#include <TGeoGlobalMagField.h>
 #include <TH1F.h> 
-#include <TF1.h> 
 #include <TRandom.h> 
+#include <TTree.h>
+#include <TVirtualMC.h> 
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
 #include "AliMagF.h"
 #include "AliPHOS.h"
-#include "AliPHOSGetter.h"
+#include "AliPHOSLoader.h"
 #include "AliRun.h"
+#include "AliRawReader.h"
 #include "AliPHOSDigitizer.h"
 #include "AliPHOSSDigitizer.h"
 #include "AliPHOSDigit.h"
@@ -109,11 +95,17 @@ class TFile;
 #include "AliPHOSCalibData.h"
 #include "AliPHOSPulseGenerator.h"
 #include "AliDAQ.h"
+#include "AliPHOSRawFitterv0.h"
+#include "AliPHOSCalibData.h"
+#include "AliPHOSRawDigiProducer.h"
+#include "AliPHOSQAChecker.h"
+#include "AliPHOSRecoParam.h"
+#include "AliPHOSSimParam.h"
 
 ClassImp(AliPHOS)
 
 //____________________________________________________________________________
-  AliPHOS:: AliPHOS() : AliDetector()
+  AliPHOS:: AliPHOS() : AliDetector(),fgCalibData(0)
 {
   // Default ctor
   fName   = "PHOS" ;
@@ -121,7 +113,8 @@ ClassImp(AliPHOS)
 }
 
 //____________________________________________________________________________
-AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title)
+AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title),
+fgCalibData(0)
 {
   //   ctor : title is used to identify the layout
 }
@@ -129,12 +122,13 @@ AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title)
 //____________________________________________________________________________
 AliPHOS::~AliPHOS() 
 {  
+  if(fgCalibData) delete fgCalibData ;
 }
 
 //____________________________________________________________________________
-AliDigitizer* AliPHOS::CreateDigitizer(AliRunDigitizer* manager) const
+AliDigitizer* AliPHOS::CreateDigitizer(AliDigitizationInput* digInput) const
 {
-  return new AliPHOSDigitizer(manager);
+  return new AliPHOSDigitizer(digInput);
 }
 
 //____________________________________________________________________________
@@ -232,7 +226,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);
 
@@ -261,7 +256,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) ;
@@ -297,9 +293,8 @@ 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() ;
+  Int_t   isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ() ;
+  Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max() ;
 
   // The scintillator of the calorimeter made of PBW04                              -> idtmed[699]
   AliMedium(0, "PHOS Xtal    $", 0, 1,
@@ -380,40 +375,11 @@ 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) ;
+}
 
-  // --- Set decent energy thresholds for gamma and electron tracking
-
-  // Tracking threshold for photons and electrons in the scintillator crystal 
-  gMC->Gstpar(idtmed[699], "CUTGAM",0.5E-4) ; 
-  gMC->Gstpar(idtmed[699], "CUTELE",1.0E-4) ;
-  // --- Generate explicitly delta rays in the titan cover ---
-  gMC->Gstpar(idtmed[704], "LOSS",3.) ;
-  gMC->Gstpar(idtmed[704], "DRAY",1.) ;
-  // --- and in aluminium parts ---
-  gMC->Gstpar(idtmed[701], "LOSS",3.) ;
-  gMC->Gstpar(idtmed[701], "DRAY",1.) ;
-  // --- and in PIN diode
-  gMC->Gstpar(idtmed[705], "LOSS",3) ;
-  gMC->Gstpar(idtmed[705], "DRAY",1) ;
-  // --- and in the passive convertor
-  gMC->Gstpar(idtmed[712], "LOSS",3) ;
-  gMC->Gstpar(idtmed[712], "DRAY",1) ;
-  // Tracking threshold for photons and electrons in the gas ArC02 
-  gMC->Gstpar(idtmed[715], "CUTGAM",1.E-5) ; 
-  gMC->Gstpar(idtmed[715], "CUTELE",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "CUTNEU",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "CUTHAD",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "CUTMUO",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "BCUTE",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "BCUTM",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "DCUTE",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "DCUTM",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "PPCUTM",1.E-5) ;
-  gMC->Gstpar(idtmed[715], "LOSS",2.) ;
-  gMC->Gstpar(idtmed[715], "DRAY",0.) ;
-  gMC->Gstpar(idtmed[715], "STRA",2.) ;
-
+//_____________________________________________________________________________
+void AliPHOS::Init()
+{
 }
 
 //____________________________________________________________________________
@@ -421,7 +387,11 @@ void AliPHOS::Digits2Raw()
 {
 // convert digits of the current event to raw data
   
-  AliPHOSLoader * loader = dynamic_cast<AliPHOSLoader*>(fLoader) ; 
+  if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
+    AliError("Energy digitization should be OFF if use Digits2Raw") ;
+  }
+
+  AliPHOSLoader * loader = static_cast<AliPHOSLoader*>(fLoader) ; 
 
   // get the digits
   loader->LoadDigits();
@@ -439,18 +409,36 @@ void AliPHOS::Digits2Raw()
     return;
   }
 
+  // get mapping from OCDB
+  const TObjArray* maps = AliPHOSRecoParam::GetMappings();
+  if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
+
   // some digitization constants
-  const Float_t    kThreshold = 0.001; // skip digits below 1 MeV
+  const Float_t    kThreshold = 1.; // skip digits below 1 ADC channel
   const Int_t      kAdcThreshold = 1;  // Lower ADC threshold to write to raw data
 
-  AliAltroBuffer* buffer = NULL;
   Int_t prevDDL = -1;
 
+  if(fgCalibData==0)
+    fgCalibData= new AliPHOSCalibData(-1) ;
+
   // Create a shaper pulse object
   AliPHOSPulseGenerator *pulse = new AliPHOSPulseGenerator();
 
+  //Set Time step of sample
+  pulse->SetTimeStep(TMath::Abs(fgCalibData->GetSampleTimeStep())) ;
+  
   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;
@@ -461,12 +449,15 @@ void AliPHOS::Digits2Raw()
 
   // loop over digits (assume ordered digits)
   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
-    AliPHOSDigit* digit = dynamic_cast<AliPHOSDigit *>(digits->At(iDigit)) ;
+    AliPHOSDigit* digit = static_cast<AliPHOSDigit *>(digits->At(iDigit)) ;
+
+    // Skip small energy below treshold
     if (digit->GetEnergy() < kThreshold) 
       continue;
+
     Int_t relId[4];
     geom->AbsToRelNumbering(digit->GetId(), relId);
-    Int_t module = relId[0];
+    Int_t module = 5-relId[0];
  
     // Begin FIXME 
     if (relId[1] != 0) 
@@ -475,88 +466,84 @@ void AliPHOS::Digits2Raw()
 
     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;
-
+    if(0<=row&&row<16 && 0<=col&&col<56) iRCU=0;
+    
     //RCU1
-    if(0<=row&&row<32 && 28<=col&&col<56) iRCU=1;
+    if(16<=row&&row<32 && 0<=col&&col<56) iRCU=1;
 
     //RCU2
-    if(32<=row&&row<64 && 0<=col&&col<28) iRCU=2;
+    if(32<=row&&row<48 && 0<=col&&col<56) iRCU=2;
 
     //RCU3
-    if(32<=row&&row<64 && 28<=col&&col<56) iRCU=3;
-
-
+    if(48<=row&&row<64 && 0<=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;
+    // here module already in PHOS online convention: 0<module<4 and inverse order.
+    Int_t iDDL = 4 * module  + iRCU;
 
     // 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);
-
-      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;
+      if (buffer[iDDL] == 0) {
+       // open new file and write dummy header
+       TString fileName = AliDAQ::DdlFileName("PHOS",iDDL);
 
+       mapping[iDDL] = (AliAltroMapping*)maps->At(iDDL);
+       buffer[iDDL]  = new AliAltroBuffer(fileName.Data(),mapping[iDDL]);
+       buffer[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
+      }
       prevDDL = iDDL;
     }
 
-    // out of time range signal (?)
-    if (digit->GetTimeR() > pulse->GetRawFormatTimeMax() ) {
-      AliInfo("Signal is out of time range.\n");
-      buffer->FillBuffer((Int_t)digit->GetEnergy());
-      buffer->FillBuffer(pulse->GetRawFormatTimeBins() );  // time bin
-      buffer->FillBuffer(3);          // bunch length      
-      buffer->WriteTrailer(3, relId[3], relId[2], module);  // trailer
+    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 ) {
+      AliDebug(2,"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 ;
-      Int_t   module = relId[0];
-      if ( digit->GetId() <= geom->GetNModules() *  geom->GetNCristalsInModule()) {
+      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;}
+       if(energy>eMax) {eMax=energy; modMax=relId[0]; colMax=col; rowMax=row;}
       }
       else {
        energy = 0; // CPV raw data format is now know yet
-      }        
+      }
       pulse->SetAmplitude(energy);
       pulse->SetTZero(digit->GetTimeR());
+      Double_t r =fgCalibData->GetHighLowRatioEmc(relId[0],relId[3],relId[2]) ;
+      pulse->SetHG2LGRatio(r) ;
       pulse->MakeSamples();
       pulse->GetSamples(adcValuesHigh, adcValuesLow) ; 
-      
-      buffer->WriteChannel(relId[3]-1, relId[2]-1, 0, 
+
+      buffer[iDDL]->WriteChannel(relId[3]-1, relId[2]-1, 0, 
                           pulse->GetRawFormatTimeBins(), adcValuesLow , kAdcThreshold);
-      buffer->WriteChannel(relId[3]-1, relId[2]-1, 1, 
+      buffer[iDDL]->WriteChannel(relId[3]-1, relId[2]-1, 1, 
                           pulse->GetRawFormatTimeBins(), adcValuesHigh, kAdcThreshold);
-      
     }
   }
-  
+  delete [] adcValuesLow;
+  delete [] adcValuesHigh;
+
   // 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",
@@ -573,9 +560,11 @@ void AliPHOS::Hits2SDigits()
 
   AliPHOSSDigitizer phosDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ;
   phosDigitizer.SetEventRange(0, -1) ; // do all the events
-  phosDigitizer.ExecuteTask("all") ; 
+  phosDigitizer.Digitize("all") ; 
 }
 
+
 //____________________________________________________________________________
 AliLoader* AliPHOS::MakeLoader(const char* topfoldername)
 {
@@ -591,9 +580,9 @@ void AliPHOS::SetTreeAddress()
   // Links Hits in the Tree to Hits array
   TBranch *branch;
   char branchname[20];
-  sprintf(branchname,"%s",GetName());
+  snprintf(branchname,20,"%s",GetName());
   // Branch address for hit tree
-    TTree *treeH = TreeH();
+    TTree *treeH = fLoader->TreeH();
   if (treeH) {
     branch = treeH->GetBranch(branchname);
     if (branch) 
@@ -605,3 +594,52 @@ void AliPHOS::SetTreeAddress()
   }
 }
 
+//____________________________________________________________________________          
+Bool_t AliPHOS::Raw2SDigits(AliRawReader* rawReader)    
+{       
+                
+  AliPHOSLoader * loader = static_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();     
+                
+  rawReader->Reset() ;
+
+  const TObjArray* maps = AliPHOSRecoParam::GetMappings();
+  if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
+
+  AliAltroMapping *mapping[20];
+  for(Int_t i = 0; i < 20; i++) {
+    mapping[i] = (AliAltroMapping*)maps->At(i);
+  }
+
+  AliPHOSRawFitterv0 fitter;
+
+  fitter.SubtractPedestals(AliPHOSSimParam::GetInstance()->EMCSubtractPedestals());
+  fitter.SetAmpOffset(AliPHOSSimParam::GetInstance()->GetGlobalAltroOffset());
+  fitter.SetAmpThreshold(AliPHOSSimParam::GetInstance()->GetGlobalAltroThreshold());
+
+  AliPHOSRawDigiProducer pr(rawReader,mapping);
+  pr.SetSampleQualityCut(AliPHOSSimParam::GetInstance()->GetEMCSampleQualityCut());     
+  pr.MakeDigits(sdigits,&fitter);
+                
+  Int_t bufferSize = 32000 ;    
+  // TBranch * sdigitsBranch = tree->Branch("PHOS",&sdigits,bufferSize);        
+  tree->Branch("PHOS",&sdigits,bufferSize);     
+  tree->Fill();         
+                
+  fLoader->WriteSDigits("OVERWRITE");   
+  return kTRUE;         
+  
+}