Fix from Yves, if fitting is not good, recalculate from parabola; reject bad channels...

diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx
index 985dfff2492bd84146923f30d55ded5fd7a3e8b5..ca9d37b24dc05f127a73e12055d5d472070666a5 100644 (file)
--- a/EMCAL/AliEMCALClusterizerv1.cxx
+++ b/EMCAL/AliEMCALClusterizerv1.cxx
@@ -184,8 +184,11 @@ Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId)
     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
          
        // Check if channel is bad (dead, hot ...), in this case return 0.      
+       // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
+       // for the moment keep it here but remember to do the selection at the sdigitizer level 
+       // and remove it from here
        if(fCaloPed->IsBadChannel(iSupMod,ieta,iphi)) {
-                 AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD!!!\n",iSupMod,ieta,iphi));
+                 AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD!!!",iSupMod,ieta,iphi));
                  return 0;
        }
          

diff --git a/EMCAL/AliEMCALRawUtils.cxx b/EMCAL/AliEMCALRawUtils.cxx
index d58056dce5770e260cf8282a9064e22be2585cd6..93b9607a1d017f003fde309fb41ea0221ad826d8 100644 (file)
--- a/EMCAL/AliEMCALRawUtils.cxx
+++ b/EMCAL/AliEMCALRawUtils.cxx
@@ -47,6 +47,7 @@ class AliCaloAltroMapping;
 class AliEMCALDigitizer;
 #include "AliEMCALDigit.h"
 #include "AliEMCAL.h"
+#include "AliCaloCalibPedestal.h"  
   
 ClassImp(AliEMCALRawUtils)
   
@@ -283,7 +284,7 @@ void AliEMCALRawUtils::Digits2Raw()
 }
 
 //____________________________________________________________________________
-void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
+void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, AliCaloCalibPedestal* pedbadmap)
 {
   // convert raw data of the current event to digits                                                                                     
 
@@ -339,6 +340,12 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
       caloFlag = in.GetCaloFlag();
       if (caloFlag != 0 && caloFlag != 1) continue; 
              
+      //Do not fit bad channels
+      if(pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) {
+       //printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow());
+       continue;
+      }  
+
       // There can be zero-suppression in the raw data, 
       // so set up the TGraph in advance
       for (i=0; i < GetRawFormatTimeBins(); i++) {
@@ -396,8 +403,9 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
        Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
        if ( (TMath::Abs(ampAsymm) > 0.1) ||
             (TMath::Abs(time - timeEstimate) > 2*GetRawFormatTimeBinWidth()) ) {
-         AliDebug(2,Form("Fit results ped %f amp %f tine %f not consistent with expectations ped %f max-ped %f time %d",
+         AliDebug(2,Form("Fit results ped %f amp %f time %f not consistent with expectations ped %f max-ped %f time %d",
                      ped, amp, time, pedEstimate, ampEstimate, timeEstimate));
+
          // what should do we do then? skip this channel or assign the simple estimate? 
          // for now just overwrite the fit results with the simple estimate
          amp = ampEstimate;
@@ -638,6 +646,62 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, const Int_t lastTimeB
     time = signalF->GetParameter(1) * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger?
     ped = signalF->GetParameter(4); 
 
+    //BEG YS alternative methods to calculate the amplitude
+    Double_t * ymx = gSig->GetX() ; 
+    Double_t * ymy = gSig->GetY() ; 
+    const Int_t kN = 3 ; 
+    Double_t ymMaxX[kN] = {0., 0., 0.} ; 
+    Double_t ymMaxY[kN] = {0., 0., 0.} ; 
+    Double_t ymax = 0. ; 
+      // find the maximum amplitude
+    Int_t ymiMax = 0 ;  
+    for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
+      if (ymy[ymi] > ymMaxY[0] ) {
+        ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
+        ymMaxX[0] = ymx[ymi] ;
+        ymiMax = ymi ; 
+      }
+    }
+      // find the maximum by fitting a parabola through the max and the two adjacent samples
+    if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
+      ymMaxY[1] = ymy[ymiMax+1] ;
+      ymMaxY[2] = ymy[ymiMax-1] ; 
+      ymMaxX[1] = ymx[ymiMax+1] ;
+      ymMaxX[2] = ymx[ymiMax-1] ; 
+      if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
+          //fit a parabola through the 3 points y= a+bx+x*x*x
+        Double_t sy = 0 ; 
+        Double_t sx = 0 ; 
+        Double_t sx2 = 0 ; 
+        Double_t sx3 = 0 ; 
+        Double_t sx4 = 0 ; 
+        Double_t sxy = 0 ; 
+        Double_t sx2y = 0 ; 
+       for (Int_t i = 0; i < kN ; i++) {
+          sy += ymMaxY[i] ; 
+          sx += ymMaxX[i] ;            
+          sx2 += ymMaxX[i]*ymMaxX[i] ; 
+          sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; 
+          sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; 
+          sxy += ymMaxX[i]*ymMaxY[i] ; 
+          sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ; 
+        }
+        Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2); 
+        Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
+        Double_t c  = cN / cD ; 
+        Double_t b  = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
+        Double_t a  = (sy-b*sx-c*sx2)/kN  ;
+        Double_t xmax = -b/(2*c) ; 
+        ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
+      }
+    }
+
+    Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ; 
+    if (diff > 0.1) 
+      amp = ymMaxY[0] ; 
+
+      //END YS
+
   } // ampEstimate > fNoiseThreshold
   return;
 }

diff --git a/EMCAL/AliEMCALRawUtils.h b/EMCAL/AliEMCALRawUtils.h
index 50e07ffc02402f72c84f48493af191b7245f1c28..018a71da859e6bef113ed5cc68a7fcdc031b159c 100644 (file)
--- a/EMCAL/AliEMCALRawUtils.h
+++ b/EMCAL/AliEMCALRawUtils.h
@@ -26,6 +26,7 @@ class TGraph;
 class TF1;
 class AliRawReader;
 class AliEMCALGeometry;
+class AliCaloCalibPedestal;
 
 class AliEMCALRawUtils : public TObject {
  public:
@@ -37,7 +38,7 @@ class AliEMCALRawUtils : public TObject {
   AliEMCALRawUtils& operator =(const AliEMCALRawUtils& rawUtils);
 
   void Digits2Raw();
-  void Raw2Digits(AliRawReader *reader,TClonesArray *digitsArr);
+  void Raw2Digits(AliRawReader *reader,TClonesArray *digitsArr, AliCaloCalibPedestal* pedbadmap);
 
   void AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time);
 

diff --git a/EMCAL/AliEMCALReconstructor.cxx b/EMCAL/AliEMCALReconstructor.cxx
index 07ae64864bd9c0d3941cbd04d35697a3a7cea92b..fb9111834073fe295b4305e16acd5b48f07c8742 100644 (file)
--- a/EMCAL/AliEMCALReconstructor.cxx
+++ b/EMCAL/AliEMCALReconstructor.cxx
@@ -118,6 +118,9 @@ AliEMCALReconstructor::~AliEMCALReconstructor()
 {
   // dtor
   delete fGeom;
+  delete fgRawUtils;
+  delete fgClusterizer;
+       
   AliCodeTimer::Instance()->Print();
 } 
 
@@ -181,11 +184,11 @@ void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits
   fgRawUtils->SetNoiseThreshold(GetRecParam()->GetNoiseThreshold());
   fgRawUtils->SetNPedSamples(GetRecParam()->GetNPedSamples());
 
-  fgRawUtils->Raw2Digits(rawReader,digitsArr);
+  fgRawUtils->Raw2Digits(rawReader,digitsArr,fPedestalData);
 
   digitsTree->Fill();
-  digitsArr->Delete();
-  delete digitsArr;
+  //digitsArr->Delete(); //Do not delete digits array are not created here.
+  //delete digitsArr;
 
 }
 
@@ -445,7 +448,7 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   pid->RunPID(esd);
   delete pid;
   
-  delete digits;
+  //delete digits;
   delete clusters;
   
   // printf(" ## AliEMCALReconstructor::FillESD() is ended : ncl %i -> %i ### \n ",nClusters, nClustersNew); 
@@ -472,7 +475,7 @@ void AliEMCALReconstructor::FillMisalMatrixes(AliESDEvent* esd)const{
        }
        //Note, that owner of copied marixes will be header
        char path[255] ;
-       TGeoHMatrix * m ;
+       TGeoHMatrix * m = 0x0;
        for(Int_t sm = 0; sm < 12; sm++){
                sprintf(path,"/ALIC_1/XEN1_1/SMOD_%d",sm+1) ; //In Geometry modules numbered 1,2,.,5
                if(sm >= 10) sprintf(path,"/ALIC_1/XEN1_1/SM10_%d",sm-10+1) ;
@@ -485,7 +488,6 @@ void AliEMCALReconstructor::FillMisalMatrixes(AliESDEvent* esd)const{
                        esd->SetEMCALMatrix(NULL,sm) ;
                }
        }
-       
 }