bug fix in the vertex selection

[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRawUtils.cxx

diff --git a/EMCAL/AliEMCALRawUtils.cxx b/EMCAL/AliEMCALRawUtils.cxx
index 23fb4dd45f9ca60fc43be8a3630ce417e1a1d15c..49c1734812aa5fefdb88b00503cf96a6a4195aa4 100644 (file)
--- a/EMCAL/AliEMCALRawUtils.cxx
+++ b/EMCAL/AliEMCALRawUtils.cxx
@@ -30,6 +30,7 @@
   
 #include "TF1.h"
 #include "TGraph.h"
+#include <TRandom.h>
 class TSystem;
   
 class AliLog;
@@ -51,6 +52,7 @@ class AliEMCALDigitizer;
 ClassImp(AliEMCALRawUtils)
   
 // Signal shape parameters
+Int_t    AliEMCALRawUtils::fgTimeBins = 256;           // number of time bins for EMCAL
 Double_t AliEMCALRawUtils::fgTimeBinWidth  = 100E-9 ; // each sample is 100 ns
 Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ;   // 15 time bins ~ 1.5 musec
 
@@ -198,8 +200,8 @@ void AliEMCALRawUtils::Digits2Raw()
   for (Int_t i=0; i < nDDL; i++)
     buffers[i] = 0;
 
-  Int_t adcValuesLow[fgkTimeBins];
-  Int_t adcValuesHigh[fgkTimeBins];
+  TArrayI adcValuesLow(fgTimeBins);
+  TArrayI adcValuesHigh(fgTimeBins);
 
   // loop over digits (assume ordered digits)
   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
@@ -261,11 +263,11 @@ void AliEMCALRawUtils::Digits2Raw()
       buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM);  // trailer
       // calculate the time response function
     } else {
-      Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; 
+      Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; 
       if (lowgain) 
-       buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold);
+       buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow.GetArray(), fgThreshold);
       else 
-       buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold);
+       buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh.GetArray(), fgThreshold);
     }
   }
   
@@ -299,7 +301,7 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
 
   AliCaloRawStreamV3 in(reader,"EMCAL",fMapping);
   // Select EMCAL DDL's;
-  reader->Select("EMCAL");
+  reader->Select("EMCAL",0,43);
 
   //Updated fitting routine from 2007 beam test takes into account
   //possibility of two peaks in data and selects first one for fitting
@@ -331,6 +333,11 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
   while (in.NextDDL()) {
     while (in.NextChannel()) {
       
+               //Check if the signal  is high or low gain and then do the fit, 
+               //if it  is from TRU do not fit
+               caloFlag = in.GetCaloFlag();
+               if (caloFlag != 0 && caloFlag != 1) continue; 
+               
       // There can be zero-suppression in the raw data, 
       // so set up the TGraph in advance
       for (i=0; i < GetRawFormatTimeBins(); i++) {
@@ -338,6 +345,8 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
       }
                
       Int_t maxTime = 0;
+         Int_t min = 0x3ff; // init to 10-bit max
+         Int_t max = 0; // init to 10-bit min
       int nsamples = 0;
       while (in.NextBunch()) {
        const UShort_t *sig = in.GetSignals();
@@ -354,20 +363,21 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
        nsamples += in.GetBunchLength();
        for (i = 0; i < in.GetBunchLength(); i++) {
          time = startBin--;
-         gSig->SetPoint(time, time, sig[i]) ;
+         gSig->SetPoint(time, time, (Double_t) sig[i]) ;
+         if (max < sig[i]) max= sig[i];
+         if (min > sig[i]) min = sig[i];
        }
       } // loop over bunches
     
       if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
 
       id =  fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
-      caloFlag = in.GetCaloFlag();
       lowGain = in.IsLowGain();
 
       gSig->Set(maxTime+1);
-      FitRaw(gSig, signalF, amp, time) ; 
+       if ( (max - min) > fNoiseThreshold) FitRaw(gSig, signalF, amp, time) ; 
     
-      if (caloFlag == 0 || caloFlag == 1) { // low gain or high gain 
+      //if (caloFlag == 0 || caloFlag == 1) { // low gain or high gain 
        if (amp > 0 && amp < 2000) {  //check both high and low end of
        //result, 2000 is somewhat arbitrary - not nice with magic numbers in the code..
          AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
@@ -375,7 +385,7 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
          AddDigit(digitsArr, id, lowGain, (Int_t)amp, time);
        }
        
-      }
+      //}
 
       // Reset graph
       for (Int_t index = 0; index < gSig->GetN(); index++) {
@@ -436,7 +446,6 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain
 void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_t & time) const 
 {
   // Fits the raw signal time distribution; from AliEMCALGetter 
-
   amp = time = 0. ; 
   Double_t ped = 0;
   Int_t nPed = 0;
@@ -457,6 +466,7 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_
     ped = 10; // put some small value as first guess
   }
 
+
   Int_t maxFound = 0;
   Int_t iMax = 0;
   Float_t max = -1;
@@ -470,7 +480,7 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_
 
   for (Int_t i=fNPedSamples; i < gSig->GetN(); i++) {
     Double_t ttime, signal;
-    gSig->GetPoint(i, ttime, signal) ; 
+    gSig->GetPoint(i, ttime, signal) ;
     if (!maxFound && signal > max) {
       iMax = i;
       max = signal;
@@ -595,6 +605,11 @@ const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const
     //Double_t noise = gRandom->Gaus(0.,fgFEENoise);
     //signal = sqrt(signal*signal + noise*noise);
 
+    // March 17,09 for fast fit simulations by Alexei Pavlinov.
+    // Get from PHOS analysis. In some sense it is open questions.
+    Double_t noise = gRandom->Gaus(0.,fgFEENoise);
+    signal += noise; 
+ 
     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
     if ( adcH[iTime] > fgkRawSignalOverflow ){  // larger than 10 bits 
       adcH[iTime] = fgkRawSignalOverflow ;