]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALRawUtils.cxx
Corrected compilation options
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRawUtils.cxx
index f59bdbb9cad742268ebdc5a44a28b4bed95f720e..3b70379d44ddc850c7b284c75799f9b3116517b9 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 "AliCDBManager.h"
 #include "AliEMCALRawUtils.h"
-#include "TF1.h"
-#include "TGraph.h"
-#include <TRandom.h>
 #include "AliRun.h"
 #include "AliRunLoader.h"
 #include "AliAltroBuffer.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)
-*/
-
+using std::vector;
 ClassImp(AliEMCALRawUtils)
 
 
@@ -74,97 +66,59 @@ AliEMCALRawUtils::AliEMCALRawUtils( Algo::fitAlgorithm fitAlgo) : fNoiseThreshol
                                                                  fGeom(0), 
                                                                  fOption(""),
                                                                  fRemoveBadChannels(kFALSE),
-                                                                 fFittingAlgorithm(0),  
+                                                                 fFittingAlgorithm(fitAlgo),  
                                                                  fTimeMin(-1.),
                                                                  fTimeMax(1.),
                                                                  fUseFALTRO(kTRUE),
                                                                  fRawAnalyzer(0),
                                                                  fTriggerRawDigitMaker(0x0)
-{
+{ // ctor; set up fit algo etc
   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
+  delete fRawAnalyzer;
+  delete fTriggerRawDigitMaker;
 }
 
 
-//____________________________________________________________________________
 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() ;
@@ -174,7 +128,7 @@ void AliEMCALRawUtils::Digits2Raw()
     return;
   }
   
-  static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
+  static const Int_t nDDL = 20*2; // 20 SM for EMCal + DCal hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here  
   AliAltroBuffer* buffers[nDDL];
   for (Int_t i=0; i < nDDL; i++)
     buffers[i] = 0;
@@ -183,36 +137,28 @@ void AliEMCALRawUtils::Digits2Raw()
   TArrayI adcValuesHigh( TIMEBINS );
   
   // loop over digits (assume ordered digits)
-  for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) 
-    {
-      AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
-      if(!digit)
-       {
-         AliFatal("NULL Digit");
-       }
-      else
-       {
-         if (digit->GetAmplitude() <  AliEMCALRawResponse::GetRawFormatThreshold() ) 
-           {
-             continue;
-           }
-         //get cell indices
-         Int_t nSM = 0;
-         Int_t nIphi = 0;
-         Int_t nIeta = 0;
-         Int_t iphi = 0;
-         Int_t ieta = 0;
-         Int_t nModule = 0;
-         fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
-         fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
-      
-         //Check which is the RCU, 0 or 1, of the cell.
-         Int_t iRCU = -111;
-         //RCU0
+  for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
+    AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
+    if(!digit) {
+      AliFatal("NULL Digit");
+    } else {
+      if (digit->GetAmplitude() <  AliEMCALRawResponse::GetRawFormatThreshold() ) {
+        continue;
+      }
+      //get cell indices
+      Int_t nSM = 0;
+      Int_t nIphi = 0;
+      Int_t nIeta = 0;
+      Int_t iphi = 0;
+      Int_t ieta = 0;
+      Int_t nModule = 0;
+      fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
+      fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
+    
+      //Check which is the RCU, 0 or 1, of the cell.
+      Int_t iRCU = -111;
       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
@@ -221,13 +167,12 @@ void AliEMCALRawUtils::Digits2Raw()
       
       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{
+      } else {
         if (buffers[iDDL] == 0) {      
           // open new file and write dummy header
           TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
@@ -237,7 +182,6 @@ void AliEMCALRawUtils::Digits2Raw()
           //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]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
         }
@@ -251,10 +195,10 @@ void AliEMCALRawUtils::Digits2Raw()
           buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM);  // trailer
           // calculate the time response function
         } else {
-          Bool_t lowgain = AliEMCALRawResponse::RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(), 
-                                                                  adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; 
-         
-         if (lowgain) 
+          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(),  AliEMCALRawResponse::GetRawFormatThreshold()  );
@@ -276,9 +220,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);
  
@@ -331,16 +276,15 @@ 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;}
   AliEMCALTriggerSTURawStream inSTU(reader);
   AliCaloRawStreamV3 in(reader,"EMCAL",fMapping);      
-  reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL
+  reader->Select("EMCAL",0,39); // 39 = AliEMCALGeoParams::fgkLastAltroDDL
   fTriggerRawDigitMaker->Reset();      
   fTriggerRawDigitMaker->SetIO(reader, in, inSTU, digitsTRG, trgData);
   fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later
@@ -348,12 +292,21 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
   Int_t lowGain  = 0;
   Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref.
   
+  Float_t bcTimePhaseCorr = 0; // for BC-based L1 phase correction
+  Int_t bcMod4 = (reader->GetBCID() % 4); // LHC uses 40 MHz, EMCal uses 10 MHz clock
+       
+  //AliCDBManager* man = AliCDBManager::Instance();
+  //Int_t runNumber = man->GetRun();
+       
+ Int_t runNumber = reader->GetRunNumber();
+
+  if ((runNumber >130850 ) && (bcMod4==0 || bcMod4==1)) 
+    bcTimePhaseCorr = -1e-7; // subtract 100 ns for certain BC values
+
   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()))
@@ -374,7 +327,7 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
              AliCaloFitResults res =  fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2());  
              if(res.GetAmp() >= fNoiseThreshold )
                {
-                 AddDigit(digitsArr, id, lowGain, res.GetAmp(),  res.GetTime(), res.GetChi2(),  res.GetNdf() ); 
+                 AddDigit(digitsArr, id, lowGain, res.GetAmp(),  res.GetTime()+bcTimePhaseCorr, res.GetChi2(),  res.GetNdf() ); 
                }
            }//ALTRO
          else if(fUseFALTRO)
@@ -389,7 +342,7 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
 
 
 void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr) 
-{
+{ // rm entries with LGnoHG (unphysical), out of time window, and too bad chi2
   AliEMCALDigit *digit = 0;
   Int_t n = 0;
   Int_t nDigits = digitsArr->GetEntriesFast();
@@ -418,95 +371,14 @@ 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 );
+{ // select which fitting algo should be used
+  delete fRawAnalyzer; // delete doesn't do anything if the pointer is 0x0
   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;
 }