]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALRawUtils.cxx
add protection against truncated events + coverity - Rachid
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRawUtils.cxx
index c750f993bb2f5d13eb178b03035e006bd96bdaef..1710d51322a1082fdc73ebd0086022aee7a74f0d 100644 (file)
 //    Need to add concurrent high and low-gain info in the future
 //    No pedestal is added to the raw signal.
 //*-- Author: Marco van Leeuwen (LBL)
-
+//*-- Major refactoring by Per Thomas Hille
 
 #include "AliEMCALRawUtils.h"
-#include "TF1.h"
-#include "TGraph.h"
-#include <TRandom.h>
 #include "AliRun.h"
 #include "AliRunLoader.h"
 #include "AliAltroBuffer.h"
 #include "AliCaloConstants.h"
 #include "AliCaloRawAnalyzer.h"
 #include "AliCaloRawAnalyzerFactory.h"
+#include "AliEMCALRawResponse.h"
 
 using namespace CALO;
 using namespace EMCAL;
 
-Double_t AliEMCALRawUtils::fgTimeTrigger  = 600E-9 ;   // the time of the trigger as approximately seen in the data
-Int_t    AliEMCALRawUtils::fgThreshold         = 1;
-Int_t    AliEMCALRawUtils::fgPedestalValue     = 0;  // pedestal value for digits2raw, default generate ZS data
-Double_t AliEMCALRawUtils::fgFEENoise          = 3.; // 3 ADC channels of noise (sampled)
-
 ClassImp(AliEMCALRawUtils)
 
 
@@ -71,7 +64,7 @@ AliEMCALRawUtils::AliEMCALRawUtils( Algo::fitAlgorithm fitAlgo) : fNoiseThreshol
                                                                  fGeom(0), 
                                                                  fOption(""),
                                                                  fRemoveBadChannels(kFALSE),
-                                                                 fFittingAlgorithm(0),  
+                                                                 fFittingAlgorithm(fitAlgo),  
                                                                  fTimeMin(-1.),
                                                                  fTimeMax(1.),
                                                                  fUseFALTRO(kTRUE),
@@ -79,89 +72,49 @@ AliEMCALRawUtils::AliEMCALRawUtils( Algo::fitAlgorithm fitAlgo) : fNoiseThreshol
                                                                  fTriggerRawDigitMaker(0x0)
 {
   SetFittingAlgorithm(fitAlgo);
-  // SetFittingAlgorithm(  Algo::kLMSOffline);
-
-  //Get Mapping RCU files from the AliEMCALRecParam                                 
   const TObjArray* maps = AliEMCALRecParam::GetMappings();
   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
-
-  for(Int_t i = 0; i < 4; i++) {
-    fMapping[i] = (AliAltroMapping*)maps->At(i);
-  }
-
-  //To make sure we match with the geometry in a simulation file,
-  //let's try to get it first.  If not, take the default geometry
+  for(Int_t i = 0; i < 4; i++) 
+    {
+      fMapping[i] = (AliAltroMapping*)maps->At(i);
+    }
+  
   AliRunLoader *rl = AliRunLoader::Instance();
-  if (rl && rl->GetAliRun()) {
+  if (rl && rl->GetAliRun()) 
+    {
     AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
-    if(emcal)fGeom = emcal->GetGeometry();
-    else {
+    if(emcal)
+      {
+       fGeom = emcal->GetGeometry();
+      }
+    else 
+      {
+       AliDebug(1, Form("Using default geometry in raw reco"));
+       fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
+      }
+    } 
+  else 
+    {
       AliDebug(1, Form("Using default geometry in raw reco"));
       fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
     }
-
-  } else {
-    AliDebug(1, Form("Using default geometry in raw reco"));
-    fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
-  }
-
-  if(!fGeom) AliFatal(Form("Could not get geometry!"));
-       
-  fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();
-
-}
-
-
-//____________________________________________________________________________
-AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry, Algo::fitAlgorithm fitAlgo) : //fHighLowGainFactor(16.), 
-                                                                                             //fOrder(2), 
-  //  fTau(2.35), 
-  fNoiseThreshold(3),
-  fNPedSamples(4), 
-  fGeom(pGeometry), 
-  fOption(""),
-  fRemoveBadChannels(kFALSE),fFittingAlgorithm(0),
-  fTimeMin(-1.),fTimeMax(1.),
-  fUseFALTRO(kTRUE),fRawAnalyzer(0),
-  fTriggerRawDigitMaker(0x0)
-{
-
- // Initialize with the given geometry - constructor required by HLT
-  // HLT does not use/support AliRunLoader(s) instances
-  // This is a minimum intervention solution
-  // Comment by MPloskon@lbl.gov
-  SetFittingAlgorithm(fitAlgo);
-  // SetFittingAlgorithm(  Algo::kLMSOffline);
-  //Get Mapping RCU files from the AliEMCALRecParam
-  const TObjArray* maps = AliEMCALRecParam::GetMappings();
-  if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
   
-  for(Int_t i = 0; i < 4; i++) 
-    {
-      fMapping[i] = (AliAltroMapping*)maps->At(i);
-    }
-
   if(!fGeom) AliFatal(Form("Could not get geometry!"));
-  fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();  
+  fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker();
 }
 
 
-//____________________________________________________________________________
 AliEMCALRawUtils::~AliEMCALRawUtils() 
 {
   //dtor
 }
 
 
-//____________________________________________________________________________
 void AliEMCALRawUtils::Digits2Raw()
 {
   // convert digits of the current event to raw data
   AliRunLoader *rl = AliRunLoader::Instance();
   AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
-  
-  // get the digits
   loader->LoadDigits("EMCAL");
   loader->GetEvent();
   TClonesArray* digits = loader->Digits() ;
@@ -189,7 +142,7 @@ void AliEMCALRawUtils::Digits2Raw()
        }
       else
        {
-         if (digit->GetAmplitude() < fgThreshold
+         if (digit->GetAmplitude() <  AliEMCALRawResponse::GetRawFormatThreshold() 
            {
              continue;
            }
@@ -205,37 +158,34 @@ void AliEMCALRawUtils::Digits2Raw()
       
          //Check which is the RCU, 0 or 1, of the cell.
          Int_t iRCU = -111;
-         //RCU0
-      if (0<=iphi&&iphi<8) iRCU=0; // first cable row
-      else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; 
-      //second cable row
-      //RCU1
-      else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; 
-      //second cable row
-      else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
-      
-      if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
-      
-      if (iRCU<0) 
-        Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU);
+         if (0<=iphi&&iphi<8) iRCU=0; // first cable row
+         else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; 
+         else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; 
+         //second cable row
+         else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
+         
+         if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
+         
+         if (iRCU<0) 
+           Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU);
       
-      //Which DDL?
-      Int_t iDDL = NRCUSPERMODULE*nSM + iRCU;
-      if (iDDL < 0 || iDDL >= nDDL){
-        Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
-      }
-      else{
-        if (buffers[iDDL] == 0) {      
-          // open new file and write dummy header
-          TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
-          //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C
+         //Which DDL?
+         Int_t iDDL = NRCUSPERMODULE*nSM + iRCU;
+         if (iDDL < 0 || iDDL >= nDDL){
+           Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
+         }
+         else{
+           if (buffers[iDDL] == 0) 
+             {      
+               // open new file and write dummy header
+               TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
+               //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C
           Int_t iRCUside=iRCU+(nSM%2)*2;
           //iRCU=0 and even (0) SM -> RCU0A.data   0
           //iRCU=1 and even (0) SM -> RCU1A.data   1
           //iRCU=0 and odd  (1) SM -> RCU0C.data   2
           //iRCU=1 and odd  (1) SM -> RCU1C.data   3
-          //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl;
-          buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]);
+         buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]);
           buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
         }
         
@@ -248,11 +198,13 @@ void AliEMCALRawUtils::Digits2Raw()
           buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM);  // trailer
           // calculate the time response function
         } else {
-          Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; 
-          if (lowgain) 
-            buffers[iDDL]->WriteChannel(ieta, iphi, 0, TIMEBINS, adcValuesLow.GetArray(), fgThreshold);
+          Bool_t lowgain = AliEMCALRawResponse::RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(), 
+                                                                  adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; 
+         
+         if (lowgain) 
+            buffers[iDDL]->WriteChannel(ieta, iphi, 0, TIMEBINS, adcValuesLow.GetArray(),  AliEMCALRawResponse::GetRawFormatThreshold()  );
           else 
-            buffers[iDDL]->WriteChannel(ieta,iphi, 1, TIMEBINS, adcValuesHigh.GetArray(), fgThreshold);
+            buffers[iDDL]->WriteChannel(ieta,iphi, 1, TIMEBINS, adcValuesHigh.GetArray(),  AliEMCALRawResponse::GetRawFormatThreshold()  );
         }
       }// iDDL under the limits
     }//digit exists
@@ -271,9 +223,10 @@ void AliEMCALRawUtils::Digits2Raw()
 }
 
 
-//____________________________________________________________________________ 
+
 void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf) 
 {
+  // comment
   AliEMCALDigit *digit = 0, *tmpdigit = 0;
   TIter nextdigit(digitsArr);
  
@@ -326,10 +279,9 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain
 }
 
 
-//____________________________________________________________________________
 void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap, TClonesArray *digitsTRG, AliEMCALTriggerData* trgData)
 {
-
+  //conversion of raw data to digits
   if(digitsArr) digitsArr->Clear("C"); 
   if (!digitsArr) { Error("Raw2Digits", "no digits found !");return;}
   if (!reader) {Error("Raw2Digits", "no raw reader found !");return;}
@@ -345,10 +297,8 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
   
   while (in.NextDDL()) 
     {
-      //  fprintf(fp," TP1\n");
       while (in.NextChannel()) 
        {
-         //      fprintf(fp," TP2\n");
          caloFlag = in.GetCaloFlag();
          if (caloFlag > 2) continue; // Work with ALTRO and FALTRO 
          if(caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow()))
@@ -413,92 +363,13 @@ void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr)
 }
 
 
-Double_t 
-AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
-{
-  Double_t signal = 0.;
-  Double_t tau    = par[2];
-  Double_t n      = par[3];
-  Double_t ped    = par[4];
-  Double_t xx     = ( x[0] - par[1] + tau ) / tau ;
-  if (xx <= 0) 
-    signal = ped ;  
-  else 
-    {  
-      signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ; 
-    }
-  return signal ;  
-}
-
-
-Bool_t AliEMCALRawUtils::RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH, 
-                                           Int_t * adcL, const Int_t keyErr) const 
-{
-  Bool_t lowGain = kFALSE ; 
-  TF1 signalF("signal", RawResponseFunction, 0, TIMEBINS, 5);
-  signalF.SetParameter(0, damp) ; 
-  signalF.SetParameter(1, (dtime + fgTimeTrigger)/ TIMEBINWITH) ; 
-  signalF.SetParameter(2, TAU) ; 
-  signalF.SetParameter(3, ORDER);
-  signalF.SetParameter(4, fgPedestalValue);
-       
-  Double_t signal=0.0, noise=0.0;
-  for (Int_t iTime = 0; iTime <  TIMEBINS; iTime++) {
-    signal = signalF.Eval(iTime) ;  
-    
-    if(keyErr>0) {
-      noise = gRandom->Gaus(0.,fgFEENoise);
-      signal += noise; 
-    }
-         
-    adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
-    if ( adcH[iTime] > MAXBINVALUE ){  // larger than 10 bits 
-      adcH[iTime] = MAXBINVALUE ;
-      lowGain = kTRUE ; 
-    }
-    signal /= HGLGFACTOR;
-    adcL[iTime] =  static_cast<Int_t>(signal + 0.5) ;
-    if ( adcL[iTime] > MAXBINVALUE )  // larger than 10 bits 
-      adcL[iTime] = MAXBINVALUE ;
-  }
-  return lowGain ; 
-}
-
-
-//__________________________________________________________________
-Double_t AliEMCALRawUtils::RawResponseFunctionLog(Double_t *x, Double_t *par)
-{
-  Double_t signal = 0. ;
-  Double_t tau    = par[2];
-  Double_t n      = par[3];
-  Double_t xx     = ( x[0] - par[1] + tau ) / tau ;
-
-  if (xx < 0)
-    { 
-      signal = par[0] - n*TMath::Log(TMath::Abs(xx)) + n * (1 - xx ) ;  
-    }
-  else 
-    {  
-      signal = par[0] + n*TMath::Log(xx) + n * (1 - xx ) ; 
-    }
-  return signal ;  
-}
-
-
-
-//__________________________________________________________________
 void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo)              
 {
-  // fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer(  Algo::kStandard );
   fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer( fitAlgo );
-  
-  //fRawAnalyzer = AliCaloRawAnalyzerFactory::CreateAnalyzer( kStandard );
-
   fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done, for the new methods
   fRawAnalyzer->SetOverflowCut ( OVERFLOWCUT );
   fRawAnalyzer->SetAmpCut(fNoiseThreshold);
   fRawAnalyzer->SetFitArrayCut(fNoiseThreshold);
-  //  return;
 }