]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALDigitizer.cxx
fix angle to absId mapping for extensions
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigitizer.cxx
index 86a775a009f671d46d420cf8408c62ba8d9043ab..2ed5b4cba3141ddf94772151965c8e76dce1ea5a 100644 (file)
@@ -67,6 +67,7 @@
 #include "AliEMCALCalibData.h"
 #include "AliEMCALSimParam.h"
 #include "AliEMCALRawDigit.h"
+#include "AliCaloCalibPedestal.h"
 
   namespace
   {
@@ -395,6 +396,9 @@ void AliEMCALDigitizer::Digitize(Int_t event)
     //Put Noise contribution, smear time and energy     
     Float_t timeResolution = 0;
     for(absID = 0; absID < nEMC; absID++){ // Nov 30, 2006 by PAI; was from 1 to nEMC
+        
+      if (IsDead(absID)) continue; // Don't instantiate dead digits
+        
       Float_t energy = 0 ;    
         
       // amplitude set to zero, noise will be added later
@@ -487,7 +491,7 @@ void AliEMCALDigitizer::Digitize(Int_t event)
           
        // add the noise now, no need for embedded events with real data
        if(!embed)
-         energy += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
+         energy += gRandom->Gaus(0., fPinNoise) ;
           
 
        // JLK 26-June-2008
@@ -641,7 +645,7 @@ void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, con
     fADCpedestalEC     = fCalibData->GetADCpedestal     (iSupMod,ieta,iphi);
     fADCchannelEC      = fCalibData->GetADCchannel      (iSupMod,ieta,iphi);
     fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
-    fTimeChannel       = fCalibData->GetTimeChannel     (iSupMod,ieta,iphi);
+    fTimeChannel       = fCalibData->GetTimeChannel     (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
     fTimeChannelDecal  = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);      
   }
   
@@ -653,6 +657,38 @@ void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, con
     energy =  fNADCEC ; 
 }
 
+//_____________________________________________________________________
+void AliEMCALDigitizer::Decalibrate(AliEMCALDigit *digit)
+{ 
+       // Load Geometry
+       const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
+       
+       if (!geom) {
+               AliFatal("Did not get geometry from EMCALLoader");
+               return;
+       }
+       
+       Int_t absId   = digit->GetId();
+       Int_t iSupMod = -1;
+       Int_t nModule = -1;
+       Int_t nIphi   = -1;
+       Int_t nIeta   = -1;
+       Int_t iphi    = -1;
+       Int_t ieta    = -1;
+       
+       Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
+       
+       if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
+       geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
+       
+       if (fCalibData) {
+               fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
+               float factor = 1./(fADCchannelEC/0.0162);
+               
+               *digit = *digit * factor;
+       }
+}
+
 //_____________________________________________________________________
 void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
 { 
@@ -686,7 +722,7 @@ void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const In
   if(fCalibData) {
     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
     fADCchannelEC  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
-    fTimeChannel   = fCalibData->GetTimeChannel(iSupMod,ieta,iphi);
+    fTimeChannel   = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
   }
   
   adc   = adc * fADCchannelEC - fADCpedestalEC;
@@ -736,8 +772,10 @@ void AliEMCALDigitizer::Digitize(Option_t *option)
     nEvents = fLastEvent - fFirstEvent + 1;
     Int_t ievent  = -1;
 
-    TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit",    32*96);
-    TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", 32*96);
+    AliEMCALGeometry *geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+    const Int_t nTRU = geom->GetNTotalTRU();
+    TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit",    nTRU*96);
+    TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", nTRU*96);
 
     rl->LoadSDigits("EMCAL");
     for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
@@ -779,8 +817,10 @@ void AliEMCALDigitizer::Digitize(Option_t *option)
   
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALDigitizer");
-    AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", 
-        gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
+    Float_t cputime   = gBenchmark->GetCpuTime("EMCALDigitizer");
+    Float_t avcputime = cputime;
+    if(nEvents==0) avcputime = 0 ;
+    AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
   } 
 }
 
@@ -809,25 +849,27 @@ void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* dig
   
   AliRunLoader *runLoader = AliRunLoader::Instance();
   
-  AliRun* run = runLoader->GetAliRun();
-  
   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
-  if(!emcalLoader){
-    AliFatal("Did not get the  Loader");
-  }
-  else {
-    AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(run->GetDetector("EMCAL"));
-    if(emcal){
-      AliEMCALGeometry* geom = emcal->GetGeometry();
+  if (!emcalLoader) AliFatal("Did not get the  Loader");
+  
+    const  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
       
       // build FOR from simulated digits
       // and xfer to the corresponding TRU input (mapping)
       
-      TClonesArray* digits = emcalLoader->Digits();
-      
+         TClonesArray* sdigits = emcalLoader->SDigits();
+               
+         TClonesArray *digits = (TClonesArray*)sdigits->Clone();
+               
+         AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
+               
       TIter NextDigit(digits);
       while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
       {
+        if (IsDead(digit)) continue;
+        
+             Decalibrate(digit);   
+                 
         Int_t id = digit->GetId();
         
         Int_t trgid;
@@ -851,8 +893,8 @@ void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* dig
       }
       
       if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
-      
-      Int_t    nSamples = 32;
+               
+      Int_t     nSamples = geom->GetNTotalTRU(); 
       Int_t *timeSamples = new Int_t[nSamples];
       
       NextDigit = TIter(digitsTMP);
@@ -861,7 +903,7 @@ void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* dig
         if (digit)
         {
           Int_t     id = digit->GetId();
-          Float_t time = digit->GetTime();
+          Float_t time = 50.e-9;
           
           Double_t depositedEnergy = 0.;
           for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
@@ -869,25 +911,25 @@ void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* dig
           if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
           
           // FIXME: Check digit time!
-          if (depositedEnergy)
-          {
+          if (depositedEnergy) {
+            depositedEnergy += gRandom->Gaus(0., .08);
             DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
             
-            for (Int_t j=0;j<nSamples;j++) 
-            {
-              timeSamples[j] = ((j << 12) & 0xFF000) | (timeSamples[j] & 0xFFF);
+            for (Int_t j=0;j<nSamples;j++) {
+              if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
+              timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
             }
             
             new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
+
+            if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
           }
         }
       }
       
       delete [] timeSamples;
-    }// AliEMCAL exists
-    else AliFatal("Could not get AliEMCAL");
-  }// loader exists
-  
+               
+         if (digits && digits->GetEntriesFast()) digits->Delete();
 }
 
 //____________________________________________________________________________ 
@@ -895,16 +937,17 @@ void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSam
 {
        // parameters:  
        // id: 0..95
-       const Int_t    reso = 11;      // 11-bit resolution ADC
-       const Double_t vFSR = 1;       // Full scale input voltage range
-       const Double_t dNe   = 125;     // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
+       const Int_t    reso = 12;      // 11-bit resolution ADC
+       const Double_t vFSR = 2.;      // Full scale input voltage range 2V (p-p)
+//     const Double_t dNe  = 125;     // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
+       const Double_t dNe  = 125/1.3; // F-ALTRO max V. FEE: factor 4 
        const Double_t vA   = .136e-6; // CSP output range: 0.136uV/e-
-       const Double_t rise = 40e-9;   // rise time (10-90%) of the FastOR signal before shaping
+       const Double_t rise = 50e-9;   // rise time (10-90%) of the FastOR signal before shaping
        
        const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
        
        Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
-       
+               
        TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
        signalF.SetParameter( 0,   vV ); 
        signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
@@ -914,8 +957,14 @@ void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSam
        {
                // FIXME: add noise (probably not simply Gaussian) according to DA measurements
                // probably plan an access to OCDB
-               
-               timeSamples[iTime] = int((TMath::Power(2, reso) / vFSR) * signalF.Eval(iTime * kTimeBinWidth) + 0.5);
+    Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
+    if (TMath::Abs(sig) > vFSR/2.) {
+      AliError("Signal overflow!");
+      timeSamples[iTime] = (1 << reso) - 1;
+    } else {
+      AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
+      timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
+    }
        }
 }
 
@@ -1185,3 +1234,79 @@ void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName
   else AliFatal("Loader not available");
 }
 
+//__________________________________________________________________
+Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit) 
+{
+  AliRunLoader *runLoader = AliRunLoader::Instance();
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
+  if (!emcalLoader) AliFatal("Did not get the  Loader");
+  
+  AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
+  if (!caloPed) {
+    AliWarning("Could not access pedestal data! No dead channel removal applied");
+    return kFALSE;
+  }
+  
+       // Load Geometry
+  const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
+  if (!geom) AliFatal("Did not get geometry from EMCALLoader");
+  
+  Int_t absId   = digit->GetId();
+  Int_t iSupMod = -1;
+  Int_t nModule = -1;
+  Int_t nIphi   = -1;
+  Int_t nIeta   = -1;
+  Int_t iphi    = -1;
+  Int_t ieta    = -1;
+       
+  Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
+       
+  if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
+  geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
+
+  Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
+
+  if (channelStatus == AliCaloCalibPedestal::kDead) 
+    return kTRUE;
+  else
+    return kFALSE;
+}
+
+
+//__________________________________________________________________
+Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
+{
+    AliRunLoader *runLoader = AliRunLoader::Instance();
+    AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
+    if (!emcalLoader) AliFatal("Did not get the  Loader");
+    
+    AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
+    if (!caloPed) {
+        AliWarning("Could not access pedestal data! No dead channel removal applied");
+        return kFALSE;
+    }
+    
+       // Load Geometry
+    const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
+    if (!geom) AliFatal("Did not get geometry from EMCALLoader");
+    
+    Int_t iSupMod = -1;
+    Int_t nModule = -1;
+    Int_t nIphi   = -1;
+    Int_t nIeta   = -1;
+    Int_t iphi    = -1;
+    Int_t ieta    = -1;
+       
+    Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
+       
+    if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
+    geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
+    
+    Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
+    
+    if (channelStatus == AliCaloCalibPedestal::kDead)
+        return kTRUE;
+    else
+        return kFALSE;
+}
+