]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALClusterizerv1.cxx
Fix from Yves, if fitting is not good, recalculate from parabola; reject bad channels...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv1.cxx
index ed4b4ee1a3a4de1980dbc25809a7a242d6085eb1..ca9d37b24dc05f127a73e12055d5d472070666a5 100644 (file)
@@ -76,6 +76,7 @@ class TSystem;
 #include "AliEMCALRecParam.h"
 #include "AliEMCALReconstructor.h"
 #include "AliCDBManager.h"
+#include "AliCaloCalibPedestal.h"
 
 class AliCDBStorage;
 #include "AliCDBEntry.h"
@@ -88,7 +89,7 @@ AliEMCALClusterizerv1::AliEMCALClusterizerv1()
     fGeom(0),
     fDefaultInit(kFALSE),
     fToUnfold(kFALSE),
-    fNumberOfECAClusters(0),fCalibData(0),
+    fNumberOfECAClusters(0),fCalibData(0),fCaloPed(0),
     fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
     fECAW0(0.),fTimeCut(0.),fMinECut(0.)
 {
@@ -103,7 +104,7 @@ AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
     fGeom(geometry),
     fDefaultInit(kFALSE),
     fToUnfold(kFALSE),
-    fNumberOfECAClusters(0),fCalibData(0),
+    fNumberOfECAClusters(0),fCalibData(0), fCaloPed(0),
     fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
     fECAW0(0.),fTimeCut(0.),fMinECut(0.)
 {
@@ -126,6 +127,27 @@ AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
 
 }
 
+//____________________________________________________________________________
+AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
+: AliEMCALClusterizer(),
+fGeom(geometry),
+fDefaultInit(kFALSE),
+fToUnfold(kFALSE),
+fNumberOfECAClusters(0),fCalibData(calib), fCaloPed(caloped),
+fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
+fECAW0(0.),fTimeCut(0.),fMinECut(0.)
+{
+       // ctor, geometry and calibration are initialized elsewhere.
+       
+       if (!fGeom)
+               AliFatal("Geometry not initialized.");
+               
+       if(!gMinuit)
+               gMinuit = new TMinuit(100) ;
+       
+}
+
+
 //____________________________________________________________________________
   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
 {
@@ -160,7 +182,16 @@ 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!!!",iSupMod,ieta,iphi));
+                 return 0;
+       }
+         
     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
   
@@ -187,7 +218,9 @@ void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
   //Get calibration parameters from file or digitizer default values.
   GetCalibrationParameters() ;
 
-
+  //Get dead channel map from file or digitizer default values.
+  GetCaloCalibPedestal() ;
+       
   fNumberOfECAClusters = 0;
 
   MakeClusters() ;  //only the real clusters
@@ -200,6 +233,8 @@ void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
   //Evaluate position, dispersion and other RecPoint properties for EC section                      
   for(index = 0; index < fRecPoints->GetEntries(); index++) {
       dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ;
+         //For each rec.point set the distance to the nearest bad crystal
+         dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed);
   }
 
   fRecPoints->Sort() ;
@@ -216,6 +251,8 @@ void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
 
   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
 
+  fRecPoints->Delete();
+
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALClusterizer");
     printf("Exec took %f seconds for Clusterizing", 
@@ -224,8 +261,8 @@ void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
 }
 
 //____________________________________________________________________________
-Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * RecPoint, AliEMCALDigit ** maxAt, 
-                                     Float_t* maxAtEnergy,
+Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, 
+                                     const Float_t* maxAtEnergy,
                                      Int_t nPar, Float_t * fitparameters) const
 {
   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
@@ -241,7 +278,7 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * RecPoint, AliEMCALDigit
   gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
   // To set the address of the minimization function
   TList * toMinuit = new TList();
-  toMinuit->AddAt(RecPoint,0) ;
+  toMinuit->AddAt(recPoint,0) ;
   toMinuit->AddAt(fDigitsArr,1) ;
   toMinuit->AddAt(fGeom,2) ;
 
@@ -326,7 +363,7 @@ void AliEMCALClusterizerv1::GetCalibrationParameters()
 
   //Check if calibration is stored in data base
 
-  if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
+  if(!fCalibData)
     {
       AliCDBEntry *entry = (AliCDBEntry*) 
        AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
@@ -338,13 +375,39 @@ void AliEMCALClusterizerv1::GetCalibrationParameters()
  
 }
 
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::GetCaloCalibPedestal() 
+{
+       // Set calibration parameters:
+       // if calibration database exists, they are read from database,
+       // otherwise, they are taken from digitizer.
+       //
+       // It is a user responsilibity to open CDB before reconstruction, 
+       // for example: 
+       // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
+       
+       //Check if calibration is stored in data base
+       
+       if(!fCaloPed)
+    {
+               AliCDBEntry *entry = (AliCDBEntry*) 
+               AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
+               if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
+    }
+       
+       if(!fCaloPed)
+               AliFatal("Pedestal info not found in CDB!");
+       
+}
+
+
 //____________________________________________________________________________
 void AliEMCALClusterizerv1::Init()
 {
   // Make all memory allocations which can not be done in default constructor.
   // Attach the Clusterizer task to the list of EMCAL tasks
   
-  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  AliRunLoader *rl = AliRunLoader::Instance();
   if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
     fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
   else 
@@ -365,7 +428,8 @@ void AliEMCALClusterizerv1::InitParameters()
   fTimeCut = 300e-9 ; // 300 ns time cut (to be tuned) 
 
   fCalibData               = 0 ;
-
+  fCaloPed                 = 0 ;
+       
   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
   if(!recParam) {
     AliFatal("Reconstruction parameters for EMCAL not set!");
@@ -513,30 +577,30 @@ void AliEMCALClusterizerv1::MakeUnfolding()
     Int_t index ;
     for(index = 0 ; index < numberofNotUnfolded ; index++){
 
-      AliEMCALRecPoint * RecPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
+      AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
 
       TVector3 gpos;
       Int_t absId;
-      RecPoint->GetGlobalPosition(gpos);
+      recPoint->GetGlobalPosition(gpos);
       fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
       if(absId > nModulesToUnfold)
         break ;
 
-      Int_t nMultipl = RecPoint->GetMultiplicity() ;
+      Int_t nMultipl = recPoint->GetMultiplicity() ;
       AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-      Int_t nMax = RecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
+      Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
 
       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
-        UnfoldCluster(RecPoint, nMax, maxAt, maxAtEnergy) ;
-        fRecPoints->Remove(RecPoint);
+        UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
+        fRecPoints->Remove(recPoint);
         fRecPoints->Compress() ;
         index-- ;
         fNumberOfECAClusters-- ;
         numberofNotUnfolded-- ;
       }
       else{
-        RecPoint->SetNExMax(1) ; //Only one local maximum
+        recPoint->SetNExMax(1) ; //Only one local maximum
       }
 
       delete[] maxAt ;
@@ -592,12 +656,12 @@ void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower,
   Float_t xpar=0.,zpar=0.,epar=0.  ;
 
   AliEMCALDigit * digit = 0 ;
-  Int_t * Digits = iniTower->GetDigitsList() ;
+  Int_t * digitsList = iniTower->GetDigitsList() ;
 
   Int_t iparam ;
   Int_t iDigit ;
   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-    digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(Digits[iDigit] ) ) ;
+    digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
     fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
     efit[iDigit] = 0;
 
@@ -616,7 +680,7 @@ void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower,
   // so that energy deposited in each cell is distributed between new clusters proportionally
   // to its contribution to efit
 
-  Float_t * Energies = iniTower->GetEnergiesList() ;
+  Float_t * energiesList = iniTower->GetEnergiesList() ;
   Float_t ratio ;
 
   iparam = 0 ;
@@ -626,24 +690,24 @@ void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower,
     epar = fitparameters[iparam+2] ;
     iparam += 3 ;
 
-    AliEMCALRecPoint * RecPoint = 0 ;
+    AliEMCALRecPoint * recPoint = 0 ;
 
     if(fNumberOfECAClusters >= fRecPoints->GetSize())
       fRecPoints->Expand(2*fNumberOfECAClusters) ;
 
     (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
-    RecPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
+    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
     fNumberOfECAClusters++ ;
-    RecPoint->SetNExMax((Int_t)nPar/3) ;
+    recPoint->SetNExMax((Int_t)nPar/3) ;
 
     Float_t eDigit ;
     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-      digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( Digits[iDigit] ) ) ;
+      digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
       fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
 
       ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
-      eDigit = Energies[iDigit] * ratio ;
-      RecPoint->AddDigit( *digit, eDigit ) ;
+      eDigit = energiesList[iDigit] * ratio ;
+      recPoint->AddDigit( *digit, eDigit ) ;
     }
   }
 
@@ -662,17 +726,17 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
 
   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 
-  AliEMCALRecPoint * RecPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;
+  AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;
   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
   // A bit buggy way to get an access to the geometry
   // To be revised!
   AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));
 
-  Int_t * Digits     = RecPoint->GetDigitsList() ;
+  Int_t * digitsList     = recPoint->GetDigitsList() ;
 
-  Int_t nOdigits = RecPoint->GetDigitsMultiplicity() ;
+  Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;
 
-  Float_t * Energies = RecPoint->GetEnergiesList() ;
+  Float_t * energiesList = recPoint->GetEnergiesList() ;
 
   fret = 0. ;
   Int_t iparam ;
@@ -688,7 +752,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
 
   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 
-    digit = dynamic_cast<AliEMCALDigit*>( digits->At( Digits[iDigit] ) );
+    digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
 
     Double_t xDigit=0 ;
     Double_t zDigit=0 ;
@@ -707,7 +771,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
         efit += x[iParam] * ShowerShape(dx,dz) ;
         iParam++ ;
       }
-      Double_t sum = 2. * (efit - Energies[iDigit]) / Energies[iDigit] ; // Here we assume, that sigma = sqrt(E)
+      Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
       iParam = 0 ;
       while(iParam < nPar ){
         Double_t xpar = x[iParam] ;
@@ -740,7 +804,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
     }
 
-    fret += (efit-Energies[iDigit])*(efit-Energies[iDigit])/Energies[iDigit] ;
+    fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
     // Here we assume, that sigma = sqrt(E) 
   }
 }
@@ -758,7 +822,7 @@ void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
     TString taskName(Version()) ;
     
     printf("--------------- "); 
-    printf(taskName.Data()) ; 
+    printf("%s",taskName.Data()) ; 
     printf(" "); 
     printf("Clusterizing digits: "); 
     printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut);