Add cell rejection based on energy in sorrounding cells, 4 cells in corners with...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 16 Oct 2011 17:42:44 +0000 (17:42 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 16 Oct 2011 17:42:44 +0000 (17:42 +0000)
Move parts of code in UserExec to independendt methods, to improve readability.
Order methods in source code in alphabetical order.

PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.cxx
PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.h

index 77bee17..55fec4f 100644 (file)
 
 ClassImp(AliAnalysisTaskEMCALClusterize)
 
-//________________________________________________________________________
+//______________________________________________________________________________
 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) 
-  : AliAnalysisTaskSE(name)
-  , fGeom(0),               fGeomName("EMCAL_COMPLETEV1") 
-  , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
-  , fCalibData(0),          fPedestalData(0)
-  , fOCDBpath("raw://"),    fAccessOCDB(kFALSE)
-  , fDigitsArr(0),          fClusterArr(0),             fCaloClusterArr(0)
-  , fRecParam(0),           fClusterizer(0)
-  , fUnfolder(0),           fJustUnfold(kFALSE) 
-  , fOutputAODBranch(0),    fOutputAODBranchName("newEMCALClusters")
-  , fFillAODFile(kTRUE),    fFillAODHeader(0)
-  , fFillAODCaloCells(0),   fRun(-1)
-  , fRecoUtils(0),          fConfigName("")
-  , fCellLabels(),          fCellSecondLabels(),        fCellTime()
-  , fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
-  , fSelectCell(kFALSE),    fSelectCellMinE(0.005),     fSelectCellMinFrac(0.001)
-  , fRemoveLEDEvents(kFALSE)
+: AliAnalysisTaskSE(name)
+, fEvent(0)
+, fGeom(0),               fGeomName("EMCAL_COMPLETEV1") 
+, fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
+, fCalibData(0),          fPedestalData(0)
+, fOCDBpath("raw://"),    fAccessOCDB(kFALSE)
+, fDigitsArr(0),          fClusterArr(0),             fCaloClusterArr(0)
+, fRecParam(0),           fClusterizer(0)
+, fUnfolder(0),           fJustUnfold(kFALSE) 
+, fOutputAODBranch(0),    fOutputAODBranchName("newEMCALClusters")
+, fFillAODFile(kTRUE),    fFillAODHeader(0)
+, fFillAODCaloCells(0),   fRun(-1)
+, fRecoUtils(0),          fConfigName("")
+, fCellLabels(),          fCellSecondLabels(),        fCellTime()
+, fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
+, fSelectCell(kFALSE),    fSelectCellMinE(0.005),     fSelectCellMinFrac(0.001)
+, fRemoveLEDEvents(kFALSE), fRemoveExoticEvents(kFALSE), fRemoveExoticCells(kFALSE)
+, fExoticCellFraction(0.9), fExoticCellDiffTime(1e6),    fExoticCellMinAmplitude(0.5)
+
 {
   //ctor
   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
@@ -101,25 +104,26 @@ AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
   
 }
 
-//________________________________________________________________________
+//______________________________________________________________
 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() 
-  : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
-, fGeom(0),               fGeomName("EMCAL_COMPLETEV1") 
-, fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
-, fCalibData(0),          fPedestalData(0)
-, fOCDBpath("raw://"),    fAccessOCDB(kFALSE)
-, fDigitsArr(0),          fClusterArr(0),             fCaloClusterArr(0)
-, fRecParam(0),           fClusterizer(0)
-, fUnfolder(0),           fJustUnfold(kFALSE) 
-, fOutputAODBranch(0),    fOutputAODBranchName("newEMCALClusters")
-, fFillAODFile(kTRUE),    fFillAODHeader(0)
-, fFillAODCaloCells(0),   fRun(-1)
-, fRecoUtils(0),          fConfigName("")
-, fCellLabels(),          fCellSecondLabels(),        fCellTime()
-, fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
-, fSelectCell(kFALSE),    fSelectCellMinE(0.005),     fSelectCellMinFrac(0.001)
-, fRemoveLEDEvents(kFALSE)
-
+: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
+, fEvent(0)
+, fGeom(0),                 fGeomName("EMCAL_COMPLETEV1") 
+, fGeomMatrixSet(kFALSE),   fLoadGeomMatrices(kFALSE)
+, fCalibData(0),            fPedestalData(0)
+, fOCDBpath("raw://"),      fAccessOCDB(kFALSE)
+, fDigitsArr(0),            fClusterArr(0),             fCaloClusterArr(0)
+, fRecParam(0),             fClusterizer(0)
+, fUnfolder(0),             fJustUnfold(kFALSE) 
+, fOutputAODBranch(0),      fOutputAODBranchName("newEMCALClusters")
+, fFillAODFile(kTRUE),      fFillAODHeader(0)
+, fFillAODCaloCells(0),     fRun(-1)
+, fRecoUtils(0),            fConfigName("")
+, fCellLabels(),            fCellSecondLabels(),        fCellTime()
+, fMaxEvent(1000000000),    fDoTrackMatching(kFALSE)
+, fSelectCell(kFALSE),      fSelectCellMinE(0.005),     fSelectCellMinFrac(0.001)
+, fRemoveLEDEvents(kFALSE), fRemoveExoticEvents(kFALSE), fRemoveExoticCells(kFALSE)
+, fExoticCellFraction(0.9), fExoticCellDiffTime(1e6),    fExoticCellMinAmplitude(0.5)
 {
   // Constructor
   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
@@ -137,7 +141,7 @@ AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
 }
 
 
-//________________________________________________________________________
+//_______________________________________________________________
 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
 {
   //dtor 
@@ -163,53 +167,393 @@ AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
   
 }
 
-//-------------------------------------------------------------------
-void AliAnalysisTaskEMCALClusterize::Init()
+//_________________________________________________________________________________
+Bool_t AliAnalysisTaskEMCALClusterize::AcceptCalibrateCell(const Int_t absID,
+                                                           Float_t  & amp, 
+                                                           Double_t & time, 
+                                                           AliVCaloCells* cells) 
 {
-  //Init analysis with configuration macro if available
+  // Reject cell if criteria not passed and calibrate it
   
-  if(gROOT->LoadMacro(fConfigName) >=0){
-    printf("Configure analysis with %s\n",fConfigName.Data());
-    AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
-    fGeomName         = clus->fGeomName; 
-    fLoadGeomMatrices = clus->fLoadGeomMatrices;
-    fOCDBpath         = clus->fOCDBpath;   
-    fAccessOCDB       = clus->fAccessOCDB;
-    fRecParam         = clus->fRecParam;
-    fJustUnfold       = clus->fJustUnfold;
-    fFillAODFile      = clus->fFillAODFile;
-    fRecoUtils        = clus->fRecoUtils; 
-    fConfigName       = clus->fConfigName;
-    fMaxEvent         = clus->fMaxEvent;
-    fDoTrackMatching  = clus->fDoTrackMatching;
-    fOutputAODBranchName = clus->fOutputAODBranchName;
-    for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
+  if(absID < 0 || absID >= 24*48*fGeom->GetNumberOfSuperModules()) return kFALSE;
+  
+  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
+  fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); 
+  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);      
+  
+  // Do not include bad channels found in analysis,
+  if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
+     fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)) {
+    return kFALSE;
+  }
+  //Recalibrate energy
+  amp  = cells->GetCellAmplitude(absID);
+  if(fRecoUtils->IsRecalibrationOn()){ 
+    amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
+  }
+  
+  // Do not include cells with too low energy, nor exotic cell
+  if(amp < fRecParam->GetMinECut() ) {
+    amp = 0.;
+    return kFALSE;
+  }
+  
+  // Recalibrate time
+  time = cells->GetCellTime(absID);
+  Int_t bc = InputEvent()->GetBunchCrossNumber();
+  
+  // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
+  if (time*1e9 < 1.) time = fCellTime[absID];
+  
+  fRecoUtils->RecalibrateCellTime(absID,bc,time);
     
+  //printf("AliAnalysisTaskEMCALClusterize::AcceptCalibrateCell() - Accepted Cell id %d, amp %f, t %f\n",absID,amp,time*1.e9);
+  
+  return kTRUE;
+}
+
+
+//_________________________________________________
+Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
+{
+  //Access to OCDB stuff
+  
+  fEvent = InputEvent();
+  if (!fEvent)
+  {
+    Warning("AccessOCDB","Event not available!!!");
+    return kFALSE;
   }
   
-}  
+  if (fEvent->GetRunNumber()==fRun)
+    return kTRUE;
+  fRun = fEvent->GetRunNumber();
+  
+  if(DebugLevel() > 1 )
+    printf("AliAnalysisTaksEMCALClusterize::AccessODCD() - Begin");
+  
+  //fGeom = AliEMCALGeometry::GetInstance(fGeomName);
+  
+  AliCDBManager *cdb = AliCDBManager::Instance();
+  
+  
+  if (fOCDBpath.Length()){
+    cdb->SetDefaultStorage(fOCDBpath.Data());
+    printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data());
+  }
+  
+  cdb->SetRun(fEvent->GetRunNumber());
+  
+  //
+  // EMCAL from RAW OCDB
+  if (fOCDBpath.Contains("alien:"))
+  {
+    cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
+    cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
+  }
+  
+  TString path = cdb->GetDefaultStorage()->GetBaseFolder();
+  
+  // init parameters:
+  
+  //Get calibration parameters 
+  if(!fCalibData)
+  {
+    AliCDBEntry *entry = (AliCDBEntry*) 
+    AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+    
+    if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
+  }
+  
+  if(!fCalibData)
+    AliFatal("Calibration parameters not found in CDB!");
+  
+  //Get calibration parameters 
+  if(!fPedestalData)
+  {
+    AliCDBEntry *entry = (AliCDBEntry*) 
+    AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
+    
+    if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
+  }
+  
+  if(!fPedestalData)
+    AliFatal("Dead map not found in CDB!");
+  
+  return kTRUE;
+}
 
-//-------------------------------------------------------------------
-void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
+//_______________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
 {
-  // Init geometry, create list of output clusters
+  // Get the input event, it can depend in embedded events what you want to get
+  // Also check if the quality of the event is good if not reject it
   
-  fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;  
-  if(fOutputAODBranchName.Length()!=0){
-    fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
-    fOutputAODBranch->SetName(fOutputAODBranchName);
-    AddAODBranch("TClonesArray", &fOutputAODBranch);
+  fEvent = 0x0;
+  
+  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+  Int_t eventN = Entry();
+  if(aodIH) eventN = aodIH->GetReadEntry(); 
+  
+  if (eventN > fMaxEvent) 
+    return ;
+  
+  //printf("Clusterizer --- Event %d-- \n",eventN);
+
+  //Check if input event are embedded events
+  //If so, take output event
+  if (aodIH && aodIH->GetMergeEvents()) {
+    fEvent  = AODEvent();
+    
+    if(!aodIH->GetMergeEMCALCells()) 
+      AliFatal("Events merged but not EMCAL cells, check analysis settings!");
+    
+    if(DebugLevel() > 1){
+     printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n");
+        printf("\t InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
+               InputEvent()->GetEMCALCells()->GetNumberOfCells());
+        printf("\t MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
+               aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
+        for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
+            AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
+            if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
+        }
+        printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
+               AODEvent()->GetEMCALCells()->GetNumberOfCells());
+    }
   }
   else {
-    AliFatal("fOutputAODBranchName not set\n");
+    fEvent =  InputEvent();
+    if(fFillAODCaloCells) FillAODCaloCells();   
+    if(fFillAODHeader)    FillAODHeader();
+  }
+  
+  if (!fEvent) {
+    Error("UserExec","Event not available");
+    return ;
   }
+  
+  //-------------------------------------------------------------------------------------
+  // Reject events if LED was firing, use only for LHC11a data 
+  // Reject event if triggered by exotic cell and remove exotic cells if not triggered
+  //-------------------------------------------------------------------------------------
+  
+  if(IsLEDEvent   ()) { fEvent = 0x0 ; return ; }
+  
+  if(IsExoticEvent()) { fEvent = 0x0 ; return ; }
+  
+  //Magic line to write events to AOD filem put after event rejection
+  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
+  
+}
+
+//____________________________________________________
+void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
+{ 
+  // Recluster calocells, transform them into digits, 
+  // feed the clusterizer with them and get new list of clusters
+  
+  //In case of MC, first loop on the clusters and fill MC label to array  
+  Int_t nClusters     = fEvent->GetNumberOfCaloClusters();
+  Int_t nClustersOrg  = 0;
+
+  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+  if(aodIH && aodIH->GetEventToMerge())  //Embedding
+    nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
+  
+  for (Int_t i = 0; i < nClusters; i++)
+  {
+    AliVCluster *clus = 0;
+    if(aodIH && aodIH->GetEventToMerge()) //Embedding
+      clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
+    else      
+      clus = fEvent->GetCaloCluster(i);
+    
+    if(!clus) return;
+    
+    if(clus->IsEMCAL()){   
+      Int_t label = clus->GetLabel();
+      Int_t label2 = -1 ;
+      if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
+      UShort_t * index    = clus->GetCellsAbsId() ;
+      for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
+        fCellLabels[index[icell]]       = label;
+        fCellSecondLabels[index[icell]] = label2;
+        fCellTime[icell]                = clus->GetTOF();        
+      }
+      nClustersOrg++;
+    }
+  } 
+  
+  // Transform CaloCells into Digits
+
+  Int_t    idigit =  0;
+  Int_t    id     = -1;
+  Float_t  amp    = -1; 
+  Double_t time   = -1; 
+  
+  AliVCaloCells *cells   = fEvent->GetEMCALCells();
+    
+  TTree *digitsTree = new TTree("digitstree","digitstree");
+  digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
+  
+  for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
+  {
+    
+    // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
+    id = cells->GetCellNumber(icell);
+    Bool_t accept =  AcceptCalibrateCell(id,amp,time,cells);
+    if(  accept && IsExoticCell(id,amp,time,cells)) accept = kFALSE;
+    if( !accept ){
+      fCellLabels[id]      =-1; //reset the entry in the array for next event
+      fCellSecondLabels[id]=-1; //reset the entry in the array for next event
+      fCellTime[id]        = 0.;        
+      if( DebugLevel() > 2 ) printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
+                                    id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
+      continue;
+    }
+        
+    //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
+    new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1); 
+    fCellLabels[id]      =-1; //reset the entry in the array for next event
+    
+    idigit++;
+  }
+  
+  //Fill the tree with digits
+  digitsTree->Fill();
+  
+  //-------------------------------------------------------------------------------------
+  //Do the clusterization
+  //-------------------------------------------------------------------------------------        
+  TTree *clustersTree = new TTree("clustertree","clustertree");
+  
+  fClusterizer->SetInput(digitsTree);
+  fClusterizer->SetOutput(clustersTree);
+  fClusterizer->Digits2Clusters("");
+  
+  //-------------------------------------------------------------------------------------
+  //Transform the recpoints into AliVClusters
+  //-------------------------------------------------------------------------------------
+  
+  clustersTree->SetBranchStatus("*",0); //disable all branches
+  clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
+  
+  TBranch *branch = clustersTree->GetBranch("EMCALECARP");
+  branch->SetAddress(&fClusterArr);
+  branch->GetEntry(0);
+  
+  RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
+  
+  if(!fCaloClusterArr){ 
+    printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
+    return;    
+  }
+  
+  if( DebugLevel() > 0 ){
+    
+    printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
+
+    if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
+      printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
+             fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
+             fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
+    }
+  }
+  
+  //Reset the array with second labels for this event
+  memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
+  
+  //---CLEAN UP-----
+  fClusterizer->Clear();
+  fDigitsArr  ->Clear("C");
+  fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
+  
+  clustersTree->Delete("all");
+  digitsTree  ->Delete("all");
+  
+}
+
+//_____________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
+{
+  // Take the event clusters and unfold them
+  
+  AliVCaloCells *cells   = fEvent->GetEMCALCells();
+  Double_t cellAmplitude = 0;
+  Double_t cellTime      = 0;
+  Short_t  cellNumber    = 0;
+  Int_t    nClustersOrg  = 0;
+  
+  // Fill the array with the EMCAL clusters, copy them
+  for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
+  {
+    AliVCluster *clus = fEvent->GetCaloCluster(i);
+    if(clus->IsEMCAL()){        
+      
+      //recalibrate/remove bad channels/etc if requested
+      if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
+        continue;
+      } 
+      
+      if(fRecoUtils->IsRecalibrationOn()){
+        
+        //Calibrate cluster
+        fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
+        
+        //CalibrateCells
+        for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
+        {
+          if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
+            break;
+          
+          Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
+          fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta); 
+          fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);      
+          
+          //Do not include bad channels found in analysis?
+          if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
+             fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
+            fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
+            fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
+            fCellTime[cellNumber]        = 0.;        
+            continue;
+          }
+          
+          cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
+          
+        }// cells loop            
+      }// recalibrate
+      
+      //Cast to ESD or AOD, needed to create the cluster array
+      AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
+      AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
+      if     (esdCluster){
+        fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );   
+      }//ESD
+      else if(aodCluster){
+        fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );   
+      }//AOD
+      else 
+        Warning("UserExec()"," - Wrong CaloCluster type?");
+      nClustersOrg++;
+    }
+  }
+  
+  //Do the unfolding
+  fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
+  
+  //CLEAN-UP
+  fUnfolder->Clear();
+  
 }
 
-//________________________________________________________________________
-void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() {
+//_____________________________________________________
+void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() 
+{
   // Put calo cells in standard branch  
-  AliVEvent * event = InputEvent();
-  AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
+  AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
   
   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
@@ -236,20 +580,20 @@ void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() {
   
 }
 
-//________________________________________________________________________
-void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
+//__________________________________________________
+void AliAnalysisTaskEMCALClusterize::FillAODHeader() 
+{
   //Put event header information in standard AOD branch
   
-  AliVEvent* event      = InputEvent();
-  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
-  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
+  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
+  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
   
   Double_t pos[3]   ;
   Double_t covVtx[6];
   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
   
   AliAODHeader* header = AODEvent()->GetHeader();
-  header->SetRunNumber(event->GetRunNumber());
+  header->SetRunNumber(fEvent->GetRunNumber());
   
   if(esdevent){
     TTree* tree = fInputHandler->GetTree();
@@ -260,14 +604,14 @@ void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
   }
   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
   
-  header->SetBunchCrossNumber(event->GetBunchCrossNumber());
-  header->SetOrbitNumber(event->GetOrbitNumber());
-  header->SetPeriodNumber(event->GetPeriodNumber());
-  header->SetEventType(event->GetEventType());
+  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
+  header->SetOrbitNumber(fEvent->GetOrbitNumber());
+  header->SetPeriodNumber(fEvent->GetPeriodNumber());
+  header->SetEventType(fEvent->GetEventType());
   
   //Centrality
-  if(event->GetCentrality()){
-    header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
+  if(fEvent->GetCentrality()){
+    header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
   }
   else{
     header->SetCentrality(0);
@@ -277,8 +621,8 @@ void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
-  header->SetTriggerMask(event->GetTriggerMask()); 
-  header->SetTriggerCluster(event->GetTriggerCluster());
+  header->SetTriggerMask(fEvent->GetTriggerMask()); 
+  header->SetTriggerCluster(fEvent->GetTriggerCluster());
   if(esdevent){
     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());    
     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());    
@@ -290,18 +634,18 @@ void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());    
   }
   
-  header->SetMagneticField(event->GetMagneticField());
+  header->SetMagneticField(fEvent->GetMagneticField());
   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.); 
   
-  header->SetZDCN1Energy(event->GetZDCN1Energy());
-  header->SetZDCP1Energy(event->GetZDCP1Energy());
-  header->SetZDCN2Energy(event->GetZDCN2Energy());
-  header->SetZDCP2Energy(event->GetZDCP2Energy());
-  header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
+  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
+  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
+  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
+  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
+  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
   
-  Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
+  Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
   Float_t diamcov[3];
-  event->GetDiamondCovXY(diamcov);
+  fEvent->GetDiamondCovXY(diamcov);
   header->SetDiamond(diamxy,diamcov);
   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
@@ -309,7 +653,7 @@ void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
   //
   //
   Int_t nVertices = 1 ;/* = prim. vtx*/;
-  Int_t nCaloClus = event->GetNumberOfCaloClusters();
+  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
   
   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
   
@@ -319,7 +663,7 @@ void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
   
   // Add primary vertex. The primary tracks will be defined
   // after the loops on the composite objects (V0, cascades, kinks)
-  event->GetPrimaryVertex()->GetXYZ(pos);
+  fEvent->GetPrimaryVertex()->GetXYZ(pos);
   Float_t chi = 0;
   if      (esdevent){
     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
@@ -332,128 +676,96 @@ void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
   
   AliAODVertex * primary = new(vertices[jVertices++])
   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
-  primary->SetName(event->GetPrimaryVertex()->GetName());
-  primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
+  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
+  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
   
 }
 
-//________________________________________________________________________
-void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
+//_________________________________________
+void AliAnalysisTaskEMCALClusterize::Init()
 {
-  // Main loop
-  // Called for each event
-  
-  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
-  Int_t eventN = Entry();
-  if(aodIH) eventN = aodIH->GetReadEntry(); 
+  //Init analysis with configuration macro if available
   
-  if (eventN > fMaxEvent) {
-    AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
-    return;
+  if(gROOT->LoadMacro(fConfigName) >=0){
+    printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
+    AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
+    fGeomName         = clus->fGeomName; 
+    fLoadGeomMatrices = clus->fLoadGeomMatrices;
+    fOCDBpath         = clus->fOCDBpath;   
+    fAccessOCDB       = clus->fAccessOCDB;
+    fRecParam         = clus->fRecParam;
+    fJustUnfold       = clus->fJustUnfold;
+    fFillAODFile      = clus->fFillAODFile;
+    fRecoUtils        = clus->fRecoUtils; 
+    fConfigName       = clus->fConfigName;
+    fMaxEvent         = clus->fMaxEvent;
+    fDoTrackMatching  = clus->fDoTrackMatching;
+    fOutputAODBranchName = clus->fOutputAODBranchName;
+    for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
+    
   }
   
-  //printf("Clusterizer --- Event %d-- \n",eventN);
-  
-  //Remove the contents of output list set in the previous event 
-  fOutputAODBranch->Clear("C");
-  
-  LoadBranches();
-  
-  //Init pointers, clusterizer, ocdb
-  if(fAccessOCDB) AccessOCDB();
-  InitClusterization();
-  
-  //Get the event
-  AliVEvent   * event    = 0;
-  AliESDEvent * esdevent = 0;
-  
-  //Fill output event with header
+}  
+
+//_______________________________________________________
+void AliAnalysisTaskEMCALClusterize::InitClusterization()
+{
+  //Select clusterization/unfolding algorithm and set all the needed parameters
   
-  //Check if input event are embedded events
-  //If so, take output event
-  if (aodIH && aodIH->GetMergeEvents()) {
-    //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\n");
-    event  = AODEvent();
-    
-    if(!aodIH->GetMergeEMCALCells()) 
-      AliFatal("Events merged but not EMCAL cells, check analysis settings!");
-    
-    //    printf("InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
-    //           InputEvent()->GetEMCALCells()->GetNumberOfCells());
-    //    printf("MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
-    //           aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
-    //    for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
-    //        AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
-    //        if(sigCluster->IsEMCAL()) printf("Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
-    //    }
-    //    printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
-    //           AODEvent()->GetEMCALCells()->GetNumberOfCells());
-    
-  }
-  else {
-    event =  InputEvent();
-    esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
-    if(fFillAODCaloCells) FillAODCaloCells();   
-    if(fFillAODHeader)    FillAODHeader();
+  if (fJustUnfold){
+    // init the unfolding afterburner 
+    delete fUnfolder;
+    fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
+    return;
   }
   
-  if (!event) {
-    Error("UserExec","Event not available");
-    AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
-    return;
+  //First init the clusterizer
+  delete fClusterizer;
+  if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
+    fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
+  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) 
+    fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
+  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){ 
+    fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
+    fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
+    fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
+  } else {
+    AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
   }
   
-  //-------------------------------------------------------------------------------------
-  // Reject event if large clusters with large energy
-  // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
-  // If clusterzer NxN or V2 it does not help
-  //-------------------------------------------------------------------------------------
-  if(fRemoveLEDEvents){
-    for (Int_t i = 0; i < InputEvent()->GetNumberOfCaloClusters(); i++)
-    {
-      AliVCluster *clus = InputEvent()->GetCaloCluster(i);
-      if(clus->IsEMCAL()){    
-        
-        if ((clus->E() > 500 && clus->GetNCells() > 200 ) || clus->GetNCells() > 200) {
-          Int_t absID = clus->GetCellsAbsId()[0];
-          Int_t sm = fGeom->GetSuperModuleNumber(absID);
-          printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster : E %f, ncells %d, absId(0) %d, SM %d\n",(Int_t)Entry(), clus->E(),  clus->GetNCells(),absID, sm);
-          AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
-          return;
-        }
-      }
-    }
-    
-    // Count number of cells with energy larger than 0.1 in SM3, cut on this number
-    Int_t ncellsSM3 = 0;
-    Int_t ncellsSM4 = 0;
-    for(Int_t icell = 0; icell < event->GetEMCALCells()->GetNumberOfCells(); icell++){
-      if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
-      if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;      
-    }
-    
-    TString triggerclasses = "";
-    if(esdevent) triggerclasses = esdevent             ->GetFiredTriggerClasses();
-    else         triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
-    
-    Int_t ncellcut = 21;
-    if(triggerclasses.Contains("EMC")) ncellcut = 35;
+  //Now set the parameters
+  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
+  fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
+  fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
+  fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
+  fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
+  fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
+  fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
+  fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
+  fClusterizer->SetInputCalibrated       ( kTRUE                               );
+  fClusterizer->SetJustClusters          ( kTRUE                               );  
+  //In case of unfolding after clusterization is requested, set the corresponding parameters
+  if(fRecParam->GetUnfold()){
+    Int_t i=0;
+    for (i = 0; i < 8; i++) {
+      fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
+    }//end of loop over parameters
+    for (i = 0; i < 3; i++) {
+      fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
+      fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
+    }//end of loop over parameters
     
-    if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 ){
-      printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
-      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
-      return;
-    }
+    fClusterizer->InitClusterUnfolding();
     
-  }// Remove LED events
-  
-  //Magic line to write events to AOD filem put after event rejection
-  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
+  }// to unfold
+}
 
+//_________________________________________________
+void AliAnalysisTaskEMCALClusterize::InitGeometry()
+{
+  // Set the geometry matrix, for the first event, skip the rest
+  // Also set once the run dependent calibrations
   
-  //-------------------------------------------------------------------------------------
-  //Set the geometry matrix, for the first event, skip the rest
-  //-------------------------------------------------------------------------------------
   if(!fGeomMatrixSet){
     if(fLoadGeomMatrices){
       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
@@ -466,18 +778,18 @@ void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
       }//SM loop
     }//Load matrices
     else if(!gGeoManager){
-      Info("UserExec","Get geo matrices from data");
+      printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs  
-      if(!strcmp(event->GetName(),"AliAODEvent")) {
+      if(!strcmp(fEvent->GetName(),"AliAODEvent")) {
         if(DebugLevel() > 1) 
           Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
       }//AOD
       else {   
         if(DebugLevel() > 1) 
-          Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
-        AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
+          printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
+        AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
         if(!esd) {
-          Error("UserExec","This event does not contain ESDs?");
+          Error("InitGeometry"," - This event does not contain ESDs?");
           AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
           return;
         }
@@ -493,430 +805,152 @@ void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
     //Recover time dependent corrections, put then in recalibration histograms. Do it once
     fRecoUtils->SetRunDependentCorrections(InputEvent()->GetRunNumber());
     
-  }//first event
-  
-  //Get the list of cells needed for unfolding and reclustering
-  AliVCaloCells *cells   = event->GetEMCALCells();
-  Int_t    nClustersOrg  = 0;
-  Double_t cellAmplitude = 0;
-  Double_t cellTime      = 0;
-  Short_t  cellNumber    = 0;
+  }//first event  
   
-  //-------------------------------------------
-  //---------Unfolding clusters----------------
-  //-------------------------------------------
-  if (fJustUnfold) {
-    
-    //Fill the array with the EMCAL clusters, copy them
-    for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
-    {
-      AliVCluster *clus = event->GetCaloCluster(i);
-      if(clus->IsEMCAL()){        
-        
-        //recalibrate/remove bad channels/etc if requested
-        if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
-          continue;
-        } 
-        
-        if(fRecoUtils->IsRecalibrationOn()){
-          
-          //Calibrate cluster
-          //printf("Energy before %f ",clus->E());
-          fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
-          //printf("; after %f\n",clus->E());
-          
-          //CalibrateCells
-          for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
-          {
-            if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
-              break;
-            
-            
-            Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
-            fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta); 
-            fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);    
-            
-            //Do not include bad channels found in analysis?
-            if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
-                fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
-              fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
-              fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
-              fCellTime[cellNumber]        = 0.;        
-              continue;
-            }
-            
-            cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
-            
-          }// cells loop            
-        }// recalibrate
-        
-        //Cast to ESD or AOD, needed to create the cluster array
-        AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
-        AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
-        if     (esdCluster){
-          fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );   
-        }//ESD
-        else if(aodCluster){
-          fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );   
-        }//AOD
-        else 
-          Warning("UserExec()"," - Wrong CaloCluster type?");
-        nClustersOrg++;
-      }
-    }
-    
-    //Do the unfolding
-    fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
-    
-    //CLEAN-UP
-    fUnfolder->Clear();
-    
-  } // just unfold ESD/AOD cluster
+}
+
+//___________________________________________________________________
+Bool_t AliAnalysisTaskEMCALClusterize::IsExoticCell(const Int_t absID, 
+                                                    const Float_t ecell, 
+                                                    const Float_t tcell, 
+                                                    AliVCaloCells* cells)
+{
+  //Look to cell neighbourhood and reject if it seems exotic
   
-  //-------------------------------------------
-  //---------- Recluster cells ----------------
-  //-------------------------------------------
+  if(!fRemoveExoticCells) return kFALSE;
   
-  else{
-    //-------------------------------------------------------------------------------------
-    //Transform CaloCells into Digits
-    //-------------------------------------------------------------------------------------
-    
-    //In case of MC, first loop on the clusters and fill MC label to array
-    //.....................................................................
-    
-    //    for(Int_t j = 0; j < 24*48*11; j++)  {
-    //      if(fCellLabels[j]      !=-1) printf("Not well initialized cell %d, label1 %d\n",j,fCellLabels[j]      ) ;
-    //      if(fCellSecondLabels[j]!=-1) printf("Not well initialized cell %d, label2 %d\n",j,fCellSecondLabels[j]) ;
-    //    }
-    
-    Int_t nClusters = event->GetNumberOfCaloClusters();
-    if(aodIH && aodIH->GetEventToMerge())  //Embedding
-      nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
-    for (Int_t i = 0; i < nClusters; i++)
-    {
-      AliVCluster *clus = 0;
-      if(aodIH && aodIH->GetEventToMerge()) //Embedding
-        clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
-      else      
-        clus = event->GetCaloCluster(i);
-      
-      if(!clus) {
-        printf("AliEMCALReclusterize::UserExec() - No Clusters\n");
-        AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
-        return;
-      }
-      
-      if(clus->IsEMCAL()){   
-        //printf("Cluster Signal %d, energy %f\n",clus->GetID(),clus->E());
-        Int_t label = clus->GetLabel();
-        Int_t label2 = -1 ;
-        if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
-        UShort_t * index    = clus->GetCellsAbsId() ;
-        for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
-          fCellLabels[index[icell]]       = label;
-          fCellSecondLabels[index[icell]] = label2;
-          //printf("Clusterizer in : TOF %g\n",clus->GetTOF()*1.e9);
-
-          fCellTime[icell]                = clus->GetTOF();        
-
-          //printf("1) absID %d, label[0] %d label[1] %d\n",index[icell], fCellLabels[index[icell]],fCellSecondLabels[index[icell]]);
-        }
-        nClustersOrg++;
-      }
-    } 
-    
-    // Create digits 
-    //.................
-    Int_t    idigit =  0;
-    Int_t    id     = -1;
-    Float_t  amp    = -1; 
-    Double_t time   = -1; 
-    
-    TTree *digitsTree = new TTree("digitstree","digitstree");
-    digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
-    
-    Int_t bc = InputEvent()->GetBunchCrossNumber();
-    
-    for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
-    {
-      if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
-        break;
-      
-      time = cellTime;
-      amp  = cellAmplitude;
-      id   = cellNumber;
-      
-      Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
-      fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta); 
-      fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);  
-      
-      //Do not include bad channels found in analysis?
-      if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
-         fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
-        fCellLabels[id]      =-1; //reset the entry in the array for next event
-        fCellSecondLabels[id]=-1; //reset the entry in the array for next event
-        fCellTime[id]        = 0.;        
-        //printf("Remove channel %d\n",id);
-        continue;
-      }
-      
-      //Recalibrate?
-      if(fRecoUtils->IsRecalibrationOn()){ 
-        //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
-        amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
-      }
-            
-      // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
-      if (time*1e9 < 1.) time = fCellTime[id];
-      
-      // Recalibrate time
-      fRecoUtils->RecalibrateCellTime(id,bc,time);
-      
-//      printf("Clusterizer: Id %d, Time org %e, Time new %e; Amp org %f, Amp new %f\n",
-//             id, cells->GetTime(icell),time, cells->GetAmplitude(icell),amp);
-
-      
-      //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
-      new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1); 
-      //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
-      //else                  printf("2) Digit cell %d, no label, amp %f \n",id,amp) ;
-      fCellLabels[id]      =-1; //reset the entry in the array for next event
-      
-      //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
-      idigit++;
-    }
-    
-    //Fill the tree with digits
-    digitsTree->Fill();
-    
-    //-------------------------------------------------------------------------------------
-    //Do the clusterization
-    //-------------------------------------------------------------------------------------        
-    TTree *clustersTree = new TTree("clustertree","clustertree");
-    
-    fClusterizer->SetInput(digitsTree);
-    fClusterizer->SetOutput(clustersTree);
-    fClusterizer->Digits2Clusters("");
-    
-    //-------------------------------------------------------------------------------------
-    //Transform the recpoints into AliVClusters
-    //-------------------------------------------------------------------------------------
-    
-    clustersTree->SetBranchStatus("*",0); //disable all branches
-    clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
-    
-    TBranch *branch = clustersTree->GetBranch("EMCALECARP");
-    branch->SetAddress(&fClusterArr);
-    branch->GetEntry(0);
-    
-    RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
-    
-    if(!fCaloClusterArr){ 
-      printf("AliAnalisysTaskEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
-      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
-      return;    
-    }
-    
-    if( DebugLevel() > 0 && fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
-      printf("AliAnalisysTaskEMCALClusterize::UserExec() - Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
-             fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
-             fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
-    }
-    
-    //Reset the array with second labels for this event
-    memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
-    
-    //---CLEAN UP-----
-    fClusterizer->Clear();
-    fDigitsArr  ->Clear("C");
-    fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
-    
-    clustersTree->Delete("all");
-    digitsTree  ->Delete("all");
-  }
+  if(ecell < fExoticCellMinAmplitude) return kFALSE;
   
-  //Recalculate track-matching for the new clusters, only with ESDs
-  if(fDoTrackMatching) fRecoUtils->FindMatches(event,fCaloClusterArr,fGeom);
+  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
+  fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); 
+  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);      
   
+  //Get close cells index, energy and time, not in corners
+  Int_t absID1 = fGeom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
+  Int_t absID2 = fGeom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
+  Int_t absID3 = fGeom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
+  Int_t absID4 = fGeom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
   
-  //-------------------------------------------------------------------------------------
-  //Put the new clusters in the AOD list
-  //-------------------------------------------------------------------------------------
+  Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
+  Double_t tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
+  Bool_t   accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
   
-  Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
-  //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
-  for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
-    AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
-    //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
-    
-    //Add matched track
-    if(fDoTrackMatching){
-      Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
-      if(trackIndex >= 0){
-        newCluster->AddTrackMatched(event->GetTrack(trackIndex));
-        if(DebugLevel() > 1) 
-          Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
-      }
-    }
-    
-    //In case of new bad channels, recalculate distance to bad channels
-    if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
-      //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
-      fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
-      //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
-    }
-    
-    //    if(newCluster->GetNLabels()>0){
-    //      printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() );
-    //      UShort_t * newindex    = newCluster->GetCellsAbsId() ;
-    //      for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){
-    //       printf("\t absID %d\n",newindex[icell]);
-    //     }
-    //    }
-    
-    //    printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
-    //    printf("labels: ");
-    //    for(UInt_t ilab = 0; ilab < newCluster->GetNLabels(); ilab++)
-    //      printf("Label[%d]: %d ",ilab, newCluster->GetLabelAt(ilab)); 
-    //     printf("\n ");
-    
-    newCluster->SetID(i);
-    new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
+  accept1 = AcceptCalibrateCell(absID1,ecell1,tcell1,cells); 
+  accept2 = AcceptCalibrateCell(absID2,ecell2,tcell2,cells); 
+  accept3 = AcceptCalibrateCell(absID3,ecell3,tcell3,cells); 
+  accept4 = AcceptCalibrateCell(absID4,ecell4,tcell4,cells); 
+
+  if(DebugLevel() > 2 ){
+    printf("AliAnalysisTaksEMCALClusterize::IsExoticCell() -  Cell absID %d \n",absID);
+    printf("\t  accept1 %d, accept2 %d, accept3 %d, accept4 %d\n",
+         accept1,accept2,accept3,accept4);
+    printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n",
+           absID,absID1,absID2,absID3,absID4);
+    printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n",
+           ecell,ecell1,ecell2,ecell3,ecell4);
+    printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n",
+           tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9,
+           TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9);
   }
   
-  //if(fOutputAODBranchName.Length()!=0)
-  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
+  if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
+  if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
+  if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
+  if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
   
-  //---CLEAN UP-----
-  fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
-}      
-
-//_____________________________________________________________________
-Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
-{
-  //Access to OCDB stuff
   
-  AliVEvent * event = InputEvent();
-  if (!event)
-  {
-    Warning("AccessOCDB","Event not available!!!");
-    return kFALSE;
+  Float_t eCross = ecell1+ecell2+ecell3+ecell4;
+  if(DebugLevel() > 2 ){  
+    printf("\t eCell %f, eCross %f, 1-eCross/eCell %f\n",ecell,eCross,1-eCross/ecell);
   }
   
-  if (event->GetRunNumber()==fRun)
+  if(1-eCross/ecell > fExoticCellFraction) {
+    //cells->SetCell(icell,absID,0,1e6);
+    if(DebugLevel() > 2 ) printf("AliAnalysisTaskEMCALClusterize::IsExoticCell() - EXOTIC CELL REMOVED id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",absID,ecell,eCross,1-eCross/ecell);
     return kTRUE;
-  fRun = event->GetRunNumber();
+  }
   
-  if(DebugLevel() > 1 )
-    Info("AccessODCD()"," Begin");
+  return kFALSE;
   
-  //fGeom = AliEMCALGeometry::GetInstance(fGeomName);
+}
+
+//____________________________________________________
+Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
+{
   
-  AliCDBManager *cdb = AliCDBManager::Instance();
+  // Check if event is exotic, get an exotic cell and compare with triggered patch
+  // If there is a match remove event ... to be completed, filled with something provisional
   
+  if(!fRemoveExoticEvents) return kFALSE;
   
-  if (fOCDBpath.Length()){
-    cdb->SetDefaultStorage(fOCDBpath.Data());
-    Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
+  // Loop on cells
+  AliVCaloCells * cells = fEvent->GetEMCALCells();
+  Float_t totCellE = 0;
+  for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
+    
+    Float_t  ecell = 0 ;
+    Double_t tcell = 0 ;
+    
+    Int_t absID   = cells->GetCellNumber(icell);
+    Bool_t accept = AcceptCalibrateCell(absID,ecell,tcell,cells);
+    if(accept && !IsExoticCell(absID, ecell,tcell,cells)) totCellE += ecell;
   }
   
-  cdb->SetRun(event->GetRunNumber());
+  //  TString triggerclasses = "";
+  //  if(esdevent) triggerclasses = esdevent             ->GetFiredTriggerClasses();
+  //  else         triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
+  //    //  
+  //    printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
+  //    AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
+  //    return;
+  //  
   
-  //
-  // EMCAL from RAW OCDB
-  if (fOCDBpath.Contains("alien:"))
-  {
-    cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
-    cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
-  }
+  //printf("TotE cell %f\n",totCellE);
+  if(totCellE < 1) return kTRUE;
   
-  TString path = cdb->GetDefaultStorage()->GetBaseFolder();
+  return kFALSE;
   
-  // init parameters:
+} 
+
+//_________________________________________________
+Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent()
+{
+  //Check if event is LED
   
-  //Get calibration parameters 
-  if(!fCalibData)
-  {
-    AliCDBEntry *entry = (AliCDBEntry*) 
-    AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+  if(!fRemoveLEDEvents) return kFALSE;
     
-    if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
+  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
+  Int_t ncellsSM3 = 0;
+  Int_t ncellsSM4 = 0;
+  AliVCaloCells * cells = fEvent->GetEMCALCells();
+  for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
+    if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
+    if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;      
   }
   
-  if(!fCalibData)
-    AliFatal("Calibration parameters not found in CDB!");
+  TString triggerclasses = "";
   
-  //Get calibration parameters 
-  if(!fPedestalData)
-  {
-    AliCDBEntry *entry = (AliCDBEntry*) 
-    AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
-    
-    if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
-  }
+  AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
+  if(esdevent) triggerclasses = esdevent              ->GetFiredTriggerClasses();
+  else         triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
   
-  if(!fPedestalData)
-    AliFatal("Dead map not found in CDB!");
-  
-  return kTRUE;
-}
-
-//________________________________________________________________________________________
-void AliAnalysisTaskEMCALClusterize::InitClusterization()
-{
-  //Select clusterization/unfolding algorithm and set all the needed parameters
+  Int_t ncellcut = 21;
+  if(triggerclasses.Contains("EMC")) ncellcut = 35;
   
-  if (fJustUnfold){
-    // init the unfolding afterburner 
-    delete fUnfolder;
-    fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
-    return;
+  if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 ){
+    printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
+    AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
+    return kTRUE;
   }
   
-  //First init the clusterizer
-  delete fClusterizer;
-  if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
-    fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
-  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) 
-    fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
-  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){ 
-    fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
-    fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
-    fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
-  } else {
-    AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
-  }
+  return kFALSE;
   
-  //Now set the parameters
-  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
-  fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
-  fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
-  fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
-  fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
-  fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
-  fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
-  fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
-  fClusterizer->SetInputCalibrated       ( kTRUE                               );
-  fClusterizer->SetJustClusters          ( kTRUE                               );  
-  //In case of unfolding after clusterization is requested, set the corresponding parameters
-  if(fRecParam->GetUnfold()){
-    Int_t i=0;
-    for (i = 0; i < 8; i++) {
-      fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
-    }//end of loop over parameters
-    for (i = 0; i < 3; i++) {
-      fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
-      fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
-    }//end of loop over parameters
-    
-    fClusterizer->InitClusterUnfolding();
+} 
 
-  }// to unfold
-}
-
-//________________________________________________________________________________________
-void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
+//______________________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, 
+                                                        TObjArray *recPoints, 
+                                                        TObjArray *clusArray)
 {
   // Restore clusters from recPoints
   // Cluster energy, global position, cells and their amplitude fractions are restored
@@ -943,13 +977,13 @@ void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
       
       absIds[ncellsTrue] = digit->GetId();
       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
-            
+      
       // In case of unfolding, remove digits with unfolded energy too low      
       if(fSelectCell){
         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)  {
           
           if(DebugLevel() > 1)  {
-            printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
+            printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
           }
           
@@ -957,14 +991,14 @@ void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
           
         } // if cuts
       }// Select cells
-            
+      
       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
       ncellsTrue++;
       clusterE  +=recPoint->GetEnergiesList()[c];
       
       // In case of embedding, fill ratio with amount of signal, 
       if (aodIH && aodIH->GetMergeEvents()) {
-       
+        
         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
@@ -1005,12 +1039,11 @@ void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
     gpos.GetXYZ(g);
     
     // create a new cluster
-    //AliAODCaloCluster *clus = new AliAODCaloCluster();
     (*clusArray)[j] = new AliAODCaloCluster() ;
     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
     j++;
     clus->SetType(AliVCluster::kEMCALClusterv1);
-    clus->SetE(clusterE);//recPoint->GetEnergy());
+    clus->SetE(clusterE);
     clus->SetPosition(g);
     clus->SetNCells(ncellsTrue);
     clus->SetCellsAbsId(absIds);
@@ -1019,7 +1052,7 @@ void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower()); 
-
+    
     if(ncells == ncellsTrue){
       Float_t elipAxis[2];
       recPoint->GetElipsAxis(elipAxis);
@@ -1037,12 +1070,6 @@ void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
       fRecoUtils->RecalculateClusterPID(clus);
       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus); 
       
-//      Float_t elipAxis[2];
-//      recPoint->GetElipsAxis(elipAxis);
-//      if(ncellsTrue > 2)
-//      printf("SS, old: l0 %f, l1 %f, D %f; New l0 %f, l1 %f, D %f\n",
-//             elipAxis[0]*elipAxis[0],elipAxis[1]*elipAxis[1], recPoint->GetDispersion(),
-//             clus->GetM02(),clus->GetM20(),clus->GetDispersion());
     }
     
     //MC
@@ -1050,11 +1077,6 @@ void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
     Int_t *parentList = recPoint->GetParents(parentMult);
     clus->SetLabel(parentList, parentMult); 
     
-    //    if(parentMult)printf("RecToESD: Labels %d",parentMult);
-    //    for (Int_t iii = 0; iii < parentMult; iii++) {
-    //      printf("\t label %d\n",parentList[iii]);
-    //    }
-    
     //Write the second major contributor to each MC cluster.
     Int_t iNewLabel ;
     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
@@ -1069,14 +1091,12 @@ void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
         }
         if (iNewLabel == 1) 
         {
-          //Int_t * newLabelArray = (Int_t*)malloc((clus->GetNLabels()+1)*sizeof(Int_t)) ;
           Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
           }
           
           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
-          //if(fCellSecondLabels[idCell]>-1)printf("\t new label %d\n",fCellSecondLabels[idCell]);
           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
           delete [] newLabelArray;
         }
@@ -1084,5 +1104,88 @@ void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
     }
     
   } // recPoints loop
-    
+  
 }
+
+//____________________________________________________________
+void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
+{
+  // Init geometry, create list of output clusters
+  
+  fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;  
+  if(fOutputAODBranchName.Length()!=0){
+    fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
+    fOutputAODBranch->SetName(fOutputAODBranchName);
+    AddAODBranch("TClonesArray", &fOutputAODBranch);
+  }
+  else {
+    AliFatal("fOutputAODBranchName not set\n");
+  }
+}
+
+//_______________________________________________________
+void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+    
+  //Remove the contents of output list set in the previous event 
+  fOutputAODBranch->Clear("C");
+  
+  LoadBranches();
+  
+  //Get the event, do some checks and settings
+  CheckAndGetEvent() ;
+  
+  if (!fEvent) {
+    if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
+    return ;
+  }
+  
+  //Init pointers, clusterizer, ocdb
+  if(fAccessOCDB) AccessOCDB();
+  
+  InitClusterization();
+  
+  InitGeometry(); // only once
+  
+  // Make clusters
+  if (fJustUnfold) ClusterUnfolding();
+  else             ClusterizeCells() ;
+  
+  //Recalculate track-matching for the new clusters
+  if(fDoTrackMatching) fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
+  
+  //Put the new clusters in the AOD list
+  
+  Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
+  for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
+    AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
+    newCluster->SetID(i);
+    new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
+
+    //Add matched track
+    if(fDoTrackMatching){
+      Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
+      if(trackIndex >= 0){
+        newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
+        if(DebugLevel() > 1) 
+          printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
+      }
+    }
+    
+    //In case of new bad channels, recalculate distance to bad channels
+    if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
+      fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
+    
+    if(DebugLevel() > 1 )    
+      printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E());
+    
+  } // cluster loop
+  
+  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
+  
+  // Clean up
+  fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
+}      
+
index b76c233..d2b5f6e 100644 (file)
@@ -36,6 +36,28 @@ class AliAnalysisTaskEMCALClusterize : public AliAnalysisTaskSE {
   virtual void   UserExec(Option_t *option);
   virtual void   Init();
   virtual void   LocalInit()                                    { Init()                       ; }
+    
+  // Event methods, settings
+  void           CheckAndGetEvent();
+  
+  Bool_t         IsExoticEvent();
+  void           SwitchOnExoticEventsRemoval()                  { fRemoveExoticEvents= kTRUE   ; }
+  void           SwitchOffExoticEventsRemoval()                 { fRemoveExoticEvents= kFALSE  ; } 
+  
+  Bool_t         IsExoticCell(const Int_t absId, const Float_t ecell, 
+                              const Float_t tcell, AliVCaloCells* cells);
+  void           SwitchOnExoticCellRemoval()                    { fRemoveExoticCells = kTRUE   ; }
+  void           SwitchOffExoticCellRemoval()                   { fRemoveExoticCells = kFALSE  ; } 
+  
+  void           SetExoticCellFractionCut(Float_t f)            { fExoticCellFraction = f      ; }
+  void           SetExoticCellDiffTimeCut(Float_t dt)           { fExoticCellDiffTime = dt     ; }
+  void           SetExoticCellMinAmplitudeCut(Float_t ma)       { fExoticCellMinAmplitude = ma ; }
+  
+  Bool_t         IsLEDEvent();
+  void           SwitchOnLEDEventsRemoval()                     { fRemoveLEDEvents   = kTRUE   ; }
+  void           SwitchOffLEDEventsRemoval()                    { fRemoveLEDEvents   = kFALSE  ; } 
+
+  Bool_t         AcceptCalibrateCell(const Int_t absId, Float_t & amp, Double_t & time, AliVCaloCells* cells) ;
   
   //OCDB
   Bool_t         AccessOCDB();
@@ -44,6 +66,7 @@ class AliAnalysisTaskEMCALClusterize : public AliAnalysisTaskSE {
   void           SetOCDBPath(const char *path)                  { fOCDBpath         = path     ; }
   
   //Geometry methods
+  void           InitGeometry();
   void           SetGeometryName(TString &name)                 { fGeomName = name             ; }
   TString        GeometryName()                          const  { return fGeomName             ; }  
   void           SwitchOnLoadOwnGeometryMatrices()              { fLoadGeomMatrices = kTRUE    ; }
@@ -61,9 +84,11 @@ class AliAnalysisTaskEMCALClusterize : public AliAnalysisTaskSE {
   void           SwitchOffFillAODCaloCells()                    { fFillAODCaloCells  = kFALSE  ; } 
   
   //Algorithms settings
-  void           JustUnfold(Bool_t yesno)                       { fJustUnfold        = yesno   ; }
   AliEMCALRecParam * GetRecParam()                       const  { return fRecParam             ; }
   void           InitClusterization();
+  void           ClusterizeCells();
+  void           ClusterUnfolding();
+  void           JustUnfold(Bool_t yesno)                       { fJustUnfold        = yesno   ; }
   
   void           SetEMCALRecoUtils(AliEMCALRecoUtils * ru)      { fRecoUtils         = ru      ; }
   AliEMCALRecoUtils* GetRecoUtils()                      const  { return fRecoUtils            ; }
@@ -78,65 +103,69 @@ class AliAnalysisTaskEMCALClusterize : public AliAnalysisTaskSE {
   void           SwitchOnCellEnergySelection()                  { fSelectCell        = kTRUE   ; }
   void           SwitchOffCellEnergySelection()                 { fSelectCell        = kFALSE  ; } 
   void           SetCellCuts(Float_t e, Float_t frac)           { fSelectCellMinE    = e       ; 
-                                                                  fSelectCellMinFrac = frac    ; }
-  void           SwitchOnLEDEventsRemoval()                     { fRemoveLEDEvents   = kTRUE   ; }
-  void           SwitchOffLEDEventsRemoval()                    { fRemoveLEDEvents   = kFALSE  ; } 
-  
+                                                                  fSelectCellMinFrac = frac    ; }  
   
  private:
     
   virtual void  RecPoints2Clusters(TClonesArray *fdigitsArr, TObjArray *fRecPoints, TObjArray *clusArray);
   
+  AliVEvent             *fEvent;                   // Event 
+  
   //Geometry  
-  AliEMCALGeometry      *fGeom;             // EMCAL geometry
-  TString                fGeomName;         // Name of geometry to use.
-  TGeoHMatrix           *fGeomMatrix[10];   // Geometry matrices with alignments
-  Bool_t                 fGeomMatrixSet;    // Set geometry matrices only once, for the first event.         
-  Bool_t                 fLoadGeomMatrices; // Matrices set from configuration, not get from geometry.root or from ESDs/AODs
+  AliEMCALGeometry      *fGeom;                    // EMCAL geometry
+  TString                fGeomName;                // Name of geometry to use.
+  TGeoHMatrix           *fGeomMatrix[10];          // Geometry matrices with alignments
+  Bool_t                 fGeomMatrixSet;           // Set geometry matrices only once, for the first event.         
+  Bool_t                 fLoadGeomMatrices;        // Matrices set from configuration, not get from geometry.root or from ESDs/AODs
 
   //OCDB
-  AliEMCALCalibData     *fCalibData;        // EMCAL calib data
-  AliCaloCalibPedestal  *fPedestalData;     // EMCAL pedestal
-  TString                fOCDBpath;         // Path with OCDB location
-  Bool_t                 fAccessOCDB;       // Need to access info from OCDB (not really)   
+  AliEMCALCalibData     *fCalibData;               // EMCAL calib data
+  AliCaloCalibPedestal  *fPedestalData;            // EMCAL pedestal
+  TString                fOCDBpath;                // Path with OCDB location
+  Bool_t                 fAccessOCDB;              // Need to access info from OCDB (not really)   
 
   //Temporal arrays
-  TClonesArray          *fDigitsArr;        //-> Digits array
-  TObjArray             *fClusterArr;       //-> Recpoints array
-  TObjArray             *fCaloClusterArr;   //-> CaloClusters array
+  TClonesArray          *fDigitsArr;               //-> Digits array
+  TObjArray             *fClusterArr;              //-> Recpoints array
+  TObjArray             *fCaloClusterArr;          //-> CaloClusters array
 
   //Clusterizers 
-  AliEMCALRecParam      *fRecParam;         // Reconstruction parameters container
-  AliEMCALClusterizer   *fClusterizer;      //! EMCAL clusterizer
-  AliEMCALAfterBurnerUF *fUnfolder;         //! Unfolding procedure
-  Bool_t                 fJustUnfold;       // Just unfold, do not recluster
+  AliEMCALRecParam      *fRecParam;                // Reconstruction parameters container
+  AliEMCALClusterizer   *fClusterizer;             //! EMCAL clusterizer
+  AliEMCALAfterBurnerUF *fUnfolder;                //! Unfolding procedure
+  Bool_t                 fJustUnfold;              // Just unfold, do not recluster
   
   //AOD
-  TClonesArray          *fOutputAODBranch;     //! AOD Branch with output clusters  
-  TString                fOutputAODBranchName; // New of output AOD branch
-  Bool_t                 fFillAODFile;         // Fill the output AOD file with the new clusters, 
-                                               // if not they will be only available for the event they were generated
-  Bool_t                 fFillAODHeader;       // Copy header to standard branch
-  Bool_t                 fFillAODCaloCells;    // Copy calocells to standard branch
-
-  Int_t                  fRun;                 //!run number
+  TClonesArray          *fOutputAODBranch;         //! AOD Branch with output clusters  
+  TString                fOutputAODBranchName;     // New of output AOD branch
+  Bool_t                 fFillAODFile;             // Fill the output AOD file with the new clusters, 
+                                                   // if not they will be only available for the event they were generated
+  Bool_t                 fFillAODHeader;           // Copy header to standard branch
+  Bool_t                 fFillAODCaloCells;        // Copy calocells to standard branch
+
+  Int_t                  fRun;                     //!run number
   
-  AliEMCALRecoUtils*     fRecoUtils;           // Access to factorized reconstruction algorithms
-  TString                fConfigName;          // Name of analysis configuration file
+  AliEMCALRecoUtils*     fRecoUtils;               // Access to factorized reconstruction algorithms
+  TString                fConfigName;              // Name of analysis configuration file
   
   Int_t                  fCellLabels[12672];       // Array with MC label to be passed to digit. 
   Int_t                  fCellSecondLabels[12672]; // Array with Second MC label to be passed to digit. 
   Double_t               fCellTime[12672];         // Array with cluster time to be passed to digit in case of AODs 
 
-  Int_t                  fMaxEvent;            // Set a maximum event
+  Int_t                  fMaxEvent;                // Set a maximum event
   
-  Bool_t                 fDoTrackMatching;     // On/Off the matching recalulation to speed up analysis in PbPb
-  Bool_t                 fSelectCell;          // Reject cells from cluster if energy is too low and recalculate position/energy and other
-  Float_t                fSelectCellMinE;      // Min energy cell threshold, after unfolding
-  Float_t                fSelectCellMinFrac;   // Min fraction of cell energy after unfolding cut
-  Bool_t                 fRemoveLEDEvents;     // Remove LED events, use only for LHC11a and if input is V1 or V1+unfolding
+  Bool_t                 fDoTrackMatching;         // On/Off the matching recalulation to speed up analysis in PbPb
+  Bool_t                 fSelectCell;              // Reject cells from cluster if energy is too low and recalculate position/energy and other
+  Float_t                fSelectCellMinE;          // Min energy cell threshold, after unfolding
+  Float_t                fSelectCellMinFrac;       // Min fraction of cell energy after unfolding cut
+  Bool_t                 fRemoveLEDEvents;         // Remove LED events, use only for LHC11a 
+  Bool_t                 fRemoveExoticEvents;      // Remove exotic events
+  Bool_t                 fRemoveExoticCells;       // Remove exotic cells
+  Float_t                fExoticCellFraction;      // Good cell if fraction < 1-ecross/ecell
+  Float_t                fExoticCellDiffTime;      // If time of candidate to exotic and close cell is too different, it must be noisy, set amp to 0
+  Float_t                fExoticCellMinAmplitude;  // Check for exotic only if amplitud is larger than this value
   
-  ClassDef(AliAnalysisTaskEMCALClusterize, 12);
+  ClassDef(AliAnalysisTaskEMCALClusterize, 14);
 
 };