]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTasks/AliAnalysisTaskCaloFilter.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGGA / CaloTasks / AliAnalysisTaskCaloFilter.cxx
index 985bf02f03b687a8ff55eb98338cfe2b3a3f12fb..e83ca504b19e8d5cfc01ab84795fc4083c4f578a 100644 (file)
@@ -18,7 +18,8 @@
 // PHOS or both, creating the corresponing AODCaloClusters
 // and AODCaloCells.
 // Keep also the AODHeader information and the vertex.
-// Keep tracks if requested
+// Keep tracks, v0s, VZERO if requested
+// Select events containing a cluster or track avobe a given threshold
 // Copy of AliAnalysisTaskESDfilter.
 // Author: Gustavo Conesa Balbastre (INFN - Frascati)
 //////////////////////////////////////////////////////////
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
 #include "AliVEventHandler.h"
+#include "AliAODHandler.h"
 #include "AliAnalysisManager.h"
 #include "AliInputEventHandler.h"
-#include "AliESDtrackCuts.h"
-#include "AliTriggerAnalysis.h"
 
 //EMCAL
 #include "AliEMCALRecoUtils.h"
@@ -52,17 +52,32 @@ ClassImp(AliAnalysisTaskCaloFilter)
 ////////////////////////////////////////////////////////////////////////
 
 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
-AliAnalysisTaskSE("CaloFilterTask"), //fCuts(0x0),
-fCaloFilter(0), fCorrect(kFALSE), 
-fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"), 
+AliAnalysisTaskSE("CaloFilterTask"), 
+fCaloFilter(0),           fEventSelection(), 
+fAcceptAllMBEvent(kFALSE),fMBTriggerMask(AliVEvent::kMB), 
+fCorrect(kFALSE), 
+fEMCALGeo(0x0),           fEMCALGeoName("EMCAL_COMPLETE12SMV1"), 
 fEMCALRecoUtils(new AliEMCALRecoUtils),
 fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
 fGeoMatrixSet(kFALSE),
-fConfigName(""),fFillAODFile(kTRUE), fFillTracks(kFALSE),
-fEnergyCut(0.), fNcellsCut (0), fVzCut(100.)
+fConfigName(""),          fFillAODFile(kTRUE), 
+fFillMCParticles(kFALSE),
+fFillTracks(kFALSE),      fFillHybridTracks(kFALSE),
+fFillAllVertices(kFALSE), fFillv0s(kFALSE),  
+fFillVZERO(kFALSE),
+fEMCALEnergyCut(0.),      fEMCALNcellsCut (0),
+fPHOSEnergyCut(0.),       fPHOSNcellsCut (0), 
+fTrackPtCut(-1),
+fVzCut(100.),             fCheckEventVertex(kTRUE),
+fEvent(0x0),              
+fESDEvent(0x0),           fAODEvent(0x0)
 {
   // Default constructor
   
+  fEventSelection[0] = kFALSE;
+  fEventSelection[1] = kFALSE;
+  fEventSelection[2] = kFALSE;
+  
   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
   //for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
   
@@ -70,17 +85,32 @@ fEnergyCut(0.), fNcellsCut (0), fVzCut(100.)
 
 //__________________________________________________
 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
-AliAnalysisTaskSE(name), //fCuts(0x0),
-fCaloFilter(0), fCorrect(kFALSE),
-fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"), 
+AliAnalysisTaskSE(name), 
+fCaloFilter(0),           fEventSelection(), 
+fAcceptAllMBEvent(kFALSE),fMBTriggerMask(AliVEvent::kMB), 
+fCorrect(kFALSE),
+fEMCALGeo(0x0),           fEMCALGeoName("EMCAL_COMPLETE12SMV1"), 
 fEMCALRecoUtils(new AliEMCALRecoUtils),
 fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
 fGeoMatrixSet(kFALSE),
-fConfigName(""),fFillAODFile(kTRUE),fFillTracks(kFALSE),
-fEnergyCut(0.), fNcellsCut (0), fVzCut(100.)
+fConfigName(""),          fFillAODFile(kTRUE),
+fFillMCParticles(kFALSE),
+fFillTracks(kFALSE),      fFillHybridTracks(kFALSE),
+fFillAllVertices(kFALSE), fFillv0s(kFALSE),
+fFillVZERO(kFALSE),
+fEMCALEnergyCut(0.),      fEMCALNcellsCut(0), 
+fPHOSEnergyCut(0.),       fPHOSNcellsCut(0), 
+fTrackPtCut(-1),
+fVzCut(100.),             fCheckEventVertex(kTRUE),
+fEvent(0x0),              
+fESDEvent(0x0),           fAODEvent(0x0)
 {
   // Constructor
   
+  fEventSelection[0] = kFALSE;
+  fEventSelection[1] = kFALSE;
+  fEventSelection[2] = kFALSE;  
+  
   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
   //for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
   
@@ -96,91 +126,132 @@ AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
   
 }
 
-//-------------------------------------------------------------------
-void AliAnalysisTaskCaloFilter::Init()
+//_____________________________________________
+Bool_t AliAnalysisTaskCaloFilter::AcceptEvent()
 {
-  //Init analysis with configuration macro if available
+  // Define conditions to accept the event to be filtered
   
-  if(gROOT->LoadMacro(fConfigName) >=0)
-  {
-    printf("Configure analysis with %s\n",fConfigName.Data());
-    
-    AliAnalysisTaskCaloFilter *filter = (AliAnalysisTaskCaloFilter*)gInterpreter->ProcessLine("ConfigCaloFilter()");
-    
-    fEMCALGeoName      = filter->fEMCALGeoName; 
-    fLoadEMCALMatrices = filter->fLoadEMCALMatrices;
-    fFillAODFile       = filter->fFillAODFile;
-    fFillTracks        = filter->fFillTracks;
-    fEMCALRecoUtils    = filter->fEMCALRecoUtils; 
-    fConfigName        = filter->fConfigName;
-    fCaloFilter        = filter->fCaloFilter;
-    fCorrect           = filter->fCorrect;
-    fEnergyCut         = filter->fEnergyCut;
-    fNcellsCut         = filter->fNcellsCut;
-    fVzCut             = filter->fVzCut;
+  if(!AcceptEventVertex()) return kFALSE;
+  
+  Bool_t eventSel = kFALSE;
+  
+  Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fMBTriggerMask);
+  
+  if     ( isMB && fAcceptAllMBEvent )                eventSel = kTRUE; // accept any MB event
+  
+  else if( fEventSelection[0] && AcceptEventEMCAL() ) eventSel = kTRUE; // accept event depending on EMCAL activity
+  
+  else if( fEventSelection[1] && AcceptEventPHOS () ) eventSel = kTRUE; // accept event depending on PHOS  activity
+  
+  else if( fEventSelection[2] && AcceptEventTrack() ) eventSel = kTRUE; // accept event depending on Track activity
     
-    for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
-  }
-} 
+  return eventSel ;
+  
+}
 
 //__________________________________________________
-void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
+Bool_t AliAnalysisTaskCaloFilter::AcceptEventEMCAL()
 {
-  // Init geometry 
-       
-  fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;  
+  // Accept event given there is a EMCAL cluster with enough energy, and not noisy, exotic
+  
+  if(fCaloFilter==kPHOS)   return kTRUE; // accept 
+  
+  if(fEMCALEnergyCut <= 0) return kTRUE; // accept
+  
+  Int_t           nCluster = InputEvent() -> GetNumberOfCaloClusters();
+  AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
+  Int_t           bc       = InputEvent() -> GetBunchCrossNumber();
+  
+  for(Int_t icalo = 0; icalo < nCluster; icalo++)
+  {
+    AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
+    
+    if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
+       fEMCALRecoUtils->IsGoodCluster(clus,fEMCALGeo,caloCell,bc))
+    {
+      
+      AliDebug(1,Form("Accept :  E %2.2f > %2.2f, nCells %d > %d",
+                      clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut));
+      
+      return kTRUE;
+    }
+    
+  }// loop
+  
+  AliDebug(1,"Reject");
+
+  //printf("Fired %s\n",((AliESDEvent*)InputEvent())->GetFiredTriggerClasses().Data());
+
+  return kFALSE;
   
 }  
 
-//____________________________________________________________
-void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
+//_________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::AcceptEventPHOS()
 {
-  // Execute analysis for current event
-  // Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
-  
-  if (fDebug > 0)  
-    printf("CaloFilter: Analysing event # %d\n", (Int_t)Entry());
+  // Accept event given there is a PHOS cluster with enough energy and not noisy/exotic
   
-  AliVEvent* event = InputEvent();
-  if(!event) 
-  {
-    printf("AliAnalysisTaskCaloFilter::UserExec - This event does not contain Input?");
-    return;
-  }
+  if(fCaloFilter==kEMCAL) return kTRUE; // accept 
   
-  // Select the event
+  if(fPHOSEnergyCut <= 0) return kTRUE; // accept
   
-  if(!AcceptEventVertex()) return;
-  if(!AcceptEventEMCAL())  return;
+  Int_t nCluster = InputEvent() -> GetNumberOfCaloClusters();
   
-  //Magic line to write events to file
-  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
+  for(Int_t icalo = 0; icalo < nCluster; icalo++)
+  {
+    AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
+    
+    if( ( clus->IsPHOS() ) && ( clus->GetNCells() > fPHOSNcellsCut ) && ( clus->E() > fPHOSEnergyCut ))
+    { 
+      
+      AliDebug(1,Form("Accept :  E %2.2f > %2.2f, nCells %d > %d",
+                      clus->E(), fPHOSEnergyCut, clus->GetNCells(), fPHOSNcellsCut));
+      
+      return kTRUE;
+    }
+    
+  }// loop
   
-  Int_t nVertices = event->GetNumberOfV0s()+event->GetNumberOfCascades();;
-  Int_t nCaloClus = event->GetNumberOfCaloClusters();
-  Int_t nTracks   = event->GetNumberOfTracks();
+  AliDebug(1,"Reject");
   
-  AODEvent()->ResetStd(nTracks, nVertices, 0, 0, 0, nCaloClus, 0, 0);
+  return kFALSE;
   
-  //
-  FillAODHeader();
+}  
+
+//__________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::AcceptEventTrack()
+{
+  // Accept event if there is a track avobe a certain pT
   
-  //
-  FillAODVertices();
+  if(fTrackPtCut <= 0) return kTRUE; // accept
   
-  //
-  FillAODTracks();
+  Double_t pTrack[3] = {0,0,0};
   
-  //
-  CorrectionsInEMCAL();
+  for (Int_t nTrack = 0; nTrack < fEvent->GetNumberOfTracks(); ++nTrack) 
+  {
+    AliVTrack *track = (AliVTrack*) fEvent->GetTrack(nTrack);
+    
+    // Select only hybrid tracks?
+    if(fAODEvent && fFillHybridTracks && !((AliAODTrack*)track)->IsHybridGlobalConstrainedGlobal()) continue;
+    
+    track->GetPxPyPz(pTrack) ;
+    
+    TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
+    
+    if(momentum.Pt() > fTrackPtCut) 
+    {
+      AliDebug(1,Form("Accept :  pT %2.2f > %2.2f",momentum.Pt(), fTrackPtCut));
+
+      return kTRUE;
+    }
+    
+  } 
   
-  //
-  FillAODCaloClusters();
+  AliDebug(1,"Reject");
   
-  //
-  FillAODCaloCells();
+  return kFALSE;
   
-}
+}  
 
 //___________________________________________________
 Bool_t AliAnalysisTaskCaloFilter::AcceptEventVertex()
@@ -190,53 +261,183 @@ Bool_t AliAnalysisTaskCaloFilter::AcceptEventVertex()
   Double_t v[3];
   InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
   
-  if(TMath::Abs(v[2]) > fVzCut) return kFALSE ;
+  if(TMath::Abs(v[2]) > fVzCut) 
+  {
+    AliDebug(1,Form("Vz Reject : vz %2.2f > %2.2f",v[2],fVzCut));
+    
+    return kFALSE ;
+  }
   
   return CheckForPrimaryVertex();
 }  
 
-//__________________________________________________
-Bool_t AliAnalysisTaskCaloFilter::AcceptEventEMCAL()
+//_______________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex()
 {
-  // Accept event given there is a cluster with enough energy
+  //Check if the vertex was well reconstructed, copy from v0Reader of conversion group
+  //It only works for ESDs
   
-  if(fCaloFilter!=kEMCAL) return kTRUE; // accept all
+  if(!fCheckEventVertex) return kTRUE;
+
+  // AODs
+  if(!fESDEvent) 
+  {
+    // Check that the vertex is not (0,0,0)
+    Double_t v[3];
+    InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
+    
+    if(TMath::Abs(v[2]) < 1e-6 && 
+       TMath::Abs(v[1]) < 1e-6 && 
+       TMath::Abs(v[0]) < 1e-6 ) 
+    {
+      AliDebug(1,"Reject v(0,0,0)");
+      
+      return kFALSE ;
+    }
+    
+    return kTRUE;
+  }
   
-  Int_t           nCluster = InputEvent() -> GetNumberOfCaloClusters();
-  AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
-  Int_t           bc       = InputEvent() -> GetBunchCrossNumber();
+  // ESDs
+  if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() > 0)
+  {
+    return kTRUE;
+  }
   
-  for(Int_t icalo = 0; icalo < nCluster; icalo++)
+  if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() < 1) 
   {
-    AliESDCaloCluster *clus = (AliESDCaloCluster*) (InputEvent()->GetCaloCluster(icalo));
-    
-    if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fNcellsCut ) && ( clus->E() > fEnergyCut ) &&
-       fEMCALRecoUtils->IsGoodCluster(clus,fEMCALGeo,caloCell,bc))
-    { 
+    // SPD vertex
+    if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() > 0) 
+    {
+      //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+      return kTRUE;
       
-      // printf("Accept event %d, E %2.2f > %2.2f, nCells %d > %d \n", 
-      //        (Int_t) Entry(),clus->E(), fEnergyCut, clus->GetNCells(), fNcellsCut);
+    }
+    if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() < 1) 
+    {
+      //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+      AliDebug(1,"Reject, GetPrimaryVertexSPD()->GetNContributors() < 1");
       
-      return kTRUE;
+      return kFALSE;
     }
-    
-  }// loop
+  }
+  
+  AliDebug(1,"Reject, GetPrimaryVertexTracks()->GetNContributors() > 1");
   
   return kFALSE;
   
-}  
+}
+
+//__________________________________________________
+void AliAnalysisTaskCaloFilter::CorrectionsInEMCAL()
+{
+  //If EMCAL, and requested, correct energy, position ...
+  
+  //Need to do this in a separate loop before filling the ESDs because of the track matching recalculations
+  
+  if(fCorrect && (fCaloFilter==kEMCAL || fCaloFilter==kBoth) ) 
+  {
+    if(!fGeoMatrixSet)
+    {
+      if(fLoadEMCALMatrices)
+      {
+        AliInfo("Load user defined EMCAL geometry matrices");
+        for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+        {
+          if(fEMCALMatrix[mod])
+          {
+            if(DebugLevel() > 1) 
+              fEMCALMatrix[mod]->Print();
+            fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;  
+          }
+          fGeoMatrixSet=kTRUE;
+        }//SM loop
+      }//Load matrices
+      else if(!gGeoManager)
+      {
+        AliInfo("Get geo matrices from data");
+        //Still not implemented in AOD, just a workaround to be able to work at least with ESDs        
+        if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) 
+        {
+          AliDebug(1,"Use ideal geometry, values geometry matrix not kept in AODs");
+        }//AOD
+        else 
+        {      
+          AliDebug(1,"Load Misaligned matrices");
+          AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
+          if(!esd) 
+          {
+            AliInfo("This event does not contain ESDs?");
+            return;
+          }
+          for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+          {
+            if(DebugLevel() > 1)
+              esd->GetEMCALMatrix(mod)->Print();
+            if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
+          } 
+          fGeoMatrixSet=kTRUE;
+        }//ESD
+      }//Load matrices from Data 
+      
+    }//first event
+    
+    //Cluster Loop
+    Int_t nCaloClus = InputEvent()->GetNumberOfCaloClusters();
+    
+    for (Int_t iClust=0; iClust<nCaloClus; ++iClust) 
+    {
+      AliVCluster * cluster = InputEvent()->GetCaloCluster(iClust);
+      
+      if(cluster->IsPHOS()) continue ;
+      
+      Float_t position[]={0,0,0};
+      
+      AliDebug(1,Form("Check cluster %d for bad channels and close to border",cluster->GetID()));
+      
+      if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;       
+      
+      AliDebug(2,Form("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f",iClust,cluster->E(),
+                      cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel()));
+      cluster->GetPosition(position);
+      AliDebug(2,Form("Filter, before  : i %d, x %f, y %f, z %f",cluster->GetID(), position[0], position[1], position[2]));
+      
+      //Recalculate distance to bad channels, if new list of bad channels provided
+      fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, InputEvent()->GetEMCALCells(), cluster);
+      
+      if(fEMCALRecoUtils->IsRecalibrationOn()) 
+      {
+        fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, InputEvent()->GetEMCALCells());
+        fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
+        fEMCALRecoUtils->RecalculateClusterPID(cluster);
+      }
+      
+      fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
+      
+      AliDebug(2,Form("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f",cluster->GetID(),cluster->E(),
+                      cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel()));
+      cluster->GetPosition(position);
+      AliDebug(1,Form("Filter, after   : i %d, x %f, y %f, z %f",cluster->GetID(), position[0], position[1], position[2]));
+      
+      cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
+      
+    }
+    
+    //Recalculate track-matching
+    fEMCALRecoUtils->FindMatches(InputEvent(),0,fEMCALGeo);
+    
+  } // corrections in EMCAL
+}
 
 //________________________________________________
 void AliAnalysisTaskCaloFilter::FillAODCaloCells()
 {  
   // Fill EMCAL/PHOS cell info
   
-  AliVEvent * event = InputEvent();
-  
   // EMCAL
-  if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && event->GetEMCALCells()) 
+  if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && fEvent->GetEMCALCells()) 
   { // protection against missing ESD information
-    AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
+    AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
     Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
     
     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
@@ -270,9 +471,9 @@ void AliAnalysisTaskCaloFilter::FillAODCaloCells()
   }
   
   // PHOS 
-  if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && event->GetPHOSCells()) 
+  if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && fEvent->GetPHOSCells()) 
   { // protection against missing ESD information
-    AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
+    AliVCaloCells &eventPHcells = *(fEvent->GetPHOSCells());
     Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
     
     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
@@ -295,19 +496,17 @@ void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
 {
   // Fill the AOD with caloclusters
   
-  AliVEvent * event = InputEvent();
-  
   // Access to the AOD container of clusters
   
   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
   Int_t jClusters=0;
   Float_t  posF[3]  ;
   
-  Int_t nCaloClus = event->GetNumberOfCaloClusters();
+  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) 
   {
     
-    AliVCluster * cluster = event->GetCaloCluster(iClust);
+    AliVCluster * cluster = fEvent->GetCaloCluster(iClust);
     
     //Check which calorimeter information we want to keep.
     
@@ -321,15 +520,13 @@ void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
     Float_t dR = cluster->GetTrackDx();
     Float_t dZ = cluster->GetTrackDz();
     
-    if(DebugLevel() > 2)
-      printf("Original residuals : dZ %f, dR %f\n ",dZ, dR);
+    AliDebug(2,Form("Original residuals : dZ %f, dR %f",dZ, dR));
     
     //--------------------------------------------------------------
     //If EMCAL and corrections done, get the new matching parameters, do not copy noisy clusters
     if(cluster->IsEMCAL() && fCorrect)
     {
-      if(DebugLevel() > 2)
-        printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
+      AliDebug(2,Form("Check cluster %d for bad channels and close to border",cluster->GetID()));
       
       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;       
       
@@ -339,11 +536,7 @@ void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
       cluster->SetTrackDistance(dR,dZ);
     }
     
-    if(DebugLevel() > 2)
-    {
-      if(cluster->IsEMCAL()) printf("EMCAL Track-Cluster Residuals : dZ %f, dR %f\n ",dZ, dR);
-      if(cluster->IsPHOS())  printf("PHOS  Track-Cluster Residuals : dZ %f, dR %f\n ",dZ, dR);
-    }
+    AliDebug(2,Form("EMCAL? %d, PHOS? %d Track-Cluster Residuals : dZ %f, dR %f",cluster->IsEMCAL(), cluster->IsPHOS(),dZ, dR));
     
     //--------------------------------------------------------------
     
@@ -374,13 +567,10 @@ void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
     caloCluster->SetTrackDistance(dR, dZ);
     
-    if(DebugLevel() > 2)
-    { 
-      printf("Filter, aod     : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
-             caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
-      caloCluster->GetPosition(posF);
-      printf("Filter, aod     : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
-    }    
+    AliDebug(2,Form("Filter, aod     : i %d, E %f, dispersion %f, m02 %f, m20 %f",caloCluster->GetID(),caloCluster->E(),
+                    caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20()));
+    caloCluster->GetPosition(posF);
+    AliDebug(2,Form("Filter, aod     : i %d, x %f, y %f, z %f",caloCluster->GetID(), posF[0], posF[1], posF[2]));
     
     //Matched tracks, just to know if there was any match, the track pointer is useless if tracks not stored
     if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990) 
@@ -394,47 +584,67 @@ void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
   
 }
 
-
-//______________________________________________
-void AliAnalysisTaskCaloFilter::FillAODHeader()
+//__________________________________________________
+void AliAnalysisTaskCaloFilter::FillAODCaloTrigger()
 {
-  // AOD header copy
+  // AOD CaloTrigger copy
   
-  AliVEvent* event = InputEvent();
-  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);  
-  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
+  if( !AODEvent() || !fAODEvent ) return;
   
-  AliAODHeader* header = AODEvent()->GetHeader();
+  AliAODCaloTrigger* triggerEM   = AODEvent()->GetCaloTrigger("EMCAL");
+  AliAODCaloTrigger* triggerPH   = AODEvent()->GetCaloTrigger("PHOS");
   
   // Copy from AODs
-  if(aodevent)
-  {
-    *header = *(aodevent->GetHeader());
-    return;
-  }
-  
-  if(!esdevent) return;
   
-  // Copy from ESDs
+  AliAODCaloTrigger* inTriggerEM = fAODEvent ->GetCaloTrigger("EMCAL");
+  AliAODCaloTrigger* inTriggerPH = fAODEvent ->GetCaloTrigger("PHOS");
   
-  header->SetRunNumber(event->GetRunNumber());
+  if(inTriggerPH && (fCaloFilter==kBoth || fCaloFilter==kPHOS))  *triggerPH = *inTriggerPH;
   
-  TTree* tree = fInputHandler->GetTree();
+  if(inTriggerEM && (fCaloFilter==kBoth || fCaloFilter==kEMCAL)) *triggerEM = *inTriggerEM;  
+}  
+
+//______________________________________________
+void AliAnalysisTaskCaloFilter::FillAODHeader()
+{
+  // AOD header copy
+  
+  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
+  if(!header)
+  {
+    AliFatal("Not a standard AOD");
+    return; // not needed but coverity complains
+  }
+  
+  // Copy from AODs
+  if(fAODEvent)
+  {
+    *header = *((AliAODHeader*)fAODEvent->GetHeader());
+    return;
+  }
+  
+  if(!fESDEvent) return;
+  
+  // Copy from ESDs
+  
+  header->SetRunNumber(fEvent->GetRunNumber());
+  
+  TTree* tree = fInputHandler->GetTree();
   if (tree) 
   {
     TFile* file = tree->GetCurrentFile();
     if (file) header->SetESDFileName(file->GetName());
   }
   
-  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())
+  if(fEvent->GetCentrality())
   {
-    header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
+    header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
   }
   else
   {
@@ -443,30 +653,45 @@ void AliAnalysisTaskCaloFilter::FillAODHeader()
   
   //Trigger  
   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
-  header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
-  header->SetTriggerMask(event->GetTriggerMask()); 
-  header->SetTriggerCluster(event->GetTriggerCluster());
-  header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());    
-  header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());    
-  header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());    
-  
-  header->SetMagneticField(event->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));
-  
-  Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
+  header->SetFiredTriggerClasses(fESDEvent->GetFiredTriggerClasses());
+  header->SetTriggerMask(fEvent->GetTriggerMask()); 
+  header->SetTriggerCluster(fEvent->GetTriggerCluster());
+  header->SetL0TriggerInputs(fESDEvent->GetHeader()->GetL0TriggerInputs());    
+  header->SetL1TriggerInputs(fESDEvent->GetHeader()->GetL1TriggerInputs());    
+  header->SetL2TriggerInputs(fESDEvent->GetHeader()->GetL2TriggerInputs());    
+  
+  header->SetMagneticField(fEvent->GetMagneticField());
+  header->SetMuonMagFieldScale(fESDEvent->GetCurrentDip()/6000.); 
+  
+  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]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
   Float_t diamcov[3];
-  event->GetDiamondCovXY(diamcov);
+  fEvent->GetDiamondCovXY(diamcov);
   header->SetDiamond(diamxy,diamcov);
-  header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
+  header->SetDiamondZ(fESDEvent->GetDiamondZ(),fESDEvent->GetSigma2DiamondZ());
   
 }
 
+
+//__________________________________________________
+void AliAnalysisTaskCaloFilter::FillAODMCParticles()
+{
+  // Copy MC particles
+  
+  if(!fFillMCParticles) return;
+  
+  TClonesArray* inMCParticles = (TClonesArray*) (fAODEvent  ->FindListObject("mcparticles"));
+  TClonesArray* ouMCParticles = (TClonesArray*) ( AODEvent()->FindListObject("mcparticles"));
+  
+  if( inMCParticles &&  ouMCParticles ) new (ouMCParticles) TClonesArray(*inMCParticles);
+    
+}  
+
 //_____________________________________________
 void AliAnalysisTaskCaloFilter::FillAODTracks()
 {
@@ -474,48 +699,179 @@ void AliAnalysisTaskCaloFilter::FillAODTracks()
   
   if(!fFillTracks) return;
   
-  AliVEvent* event = InputEvent();
-  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);  
-  //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
+  AliAODTrack* aodTrack(0x0);
   
+  Double_t pos[3]   = { 0. };      
+  Double_t covTr[21]= { 0. };
+  //Double_t pid[10]  = { 0. };  
+  Double_t p[3]     = { 0. };
+    
   // Copy from AODs
-  if(aodevent)
+  if(fAODEvent)
   {
-    TClonesArray* inTracks = aodevent  ->GetTracks();
+    //TClonesArray* inTracks = fAODEvent  ->GetTracks();
     TClonesArray* ouTracks = AODEvent()->GetTracks();
+               //new (ouTracks) TClonesArray(*inTracks);
     
-               new (ouTracks) TClonesArray(*inTracks);
+    //printf("N tracks %d\n",fAODEvent->GetNumberOfTracks());
+    Int_t nCopyTrack = 0;
+    for (Int_t nTrack = 0; nTrack < fAODEvent->GetNumberOfTracks(); ++nTrack) 
+    {
+      AliAODTrack *track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(nTrack));
+      if(!track) AliFatal("Not a standard AOD");
+      
+      // Select only hybrid tracks?
+      if(fFillHybridTracks && !track->IsHybridGlobalConstrainedGlobal()) continue;
+      
+      // Remove PID object to save space
+      //track->SetDetPID(0x0);
+      
+      //new((*ouTracks)[nCopyTrack++]) AliAODTrack(*track);
+      
+      track->GetPxPyPz(p);
+      Bool_t isDCA  = track->GetPosition(pos);
+      track->GetCovMatrix(covTr);
+      //track->GetPID(pid);
+      
+      AliAODVertex* primVertex = (AliAODVertex*) AODEvent()->GetVertices()->At(0); // primary vertex, copied previously!!!
+
+      aodTrack = new((*ouTracks)[nCopyTrack++]) AliAODTrack(
+                                                            track->GetID(),
+                                                            track->GetLabel(),
+                                                            p,
+                                                            kTRUE,
+                                                            pos,
+                                                            isDCA,
+                                                            covTr, 
+                                                            track->Charge(),
+                                                            track->GetITSClusterMap(), 
+                                                            // pid,
+                                                            primVertex,
+                                                            track->GetUsedForVtxFit(),
+                                                            track->GetUsedForPrimVtxFit(),
+                                                            (AliAODTrack::AODTrk_t) track->GetType(), 
+                                                            track->GetFilterMap());
+      
+      
+      aodTrack->SetPIDForTracking(track->GetPIDForTracking());
+      aodTrack->SetIsHybridGlobalConstrainedGlobal(track->IsHybridGlobalConstrainedGlobal());   
+      aodTrack->SetIsHybridTPCConstrainedGlobal   (track->IsHybridTPCConstrainedGlobal());    
+      aodTrack->SetIsGlobalConstrained            (track->IsGlobalConstrained());  
+      aodTrack->SetIsTPCConstrained               (track->IsTPCConstrained());    
+      
+      aodTrack->SetTPCFitMap    (track->GetTPCFitMap());
+      aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());
+      aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());
+      
+      aodTrack->SetChi2MatchTrigger(track->GetChi2MatchTrigger());
+      
+      // set the DCA values to the AOD track
+    
+      aodTrack->SetPxPyPzAtDCA(track->PxAtDCA(),track->PyAtDCA(),track->PzAtDCA());
+      aodTrack->SetXYAtDCA    (track->XAtDCA() ,track->YAtDCA());
+      
+      aodTrack->SetFlags      (track->GetFlags());
+      aodTrack->SetTPCPointsF (track->GetTPCNclsF());
+      
+      // Calo 
+      
+      if(track->IsEMCAL()) aodTrack->SetEMCALcluster(track->GetEMCALcluster());
+      if(track->IsPHOS())  aodTrack->SetPHOScluster (track->GetPHOScluster());
+      aodTrack->SetTrackPhiEtaPtOnEMCal( track->GetTrackPhiOnEMCal(),  track->GetTrackPhiOnEMCal(),  track->GetTrackPtOnEMCal() );
+          
+    } 
+    
+    //printf("Final N tracks %d\n",nCopyTrack);
+    
+    return;
+  }
+  
+}  
+
+//_________________________________________
+void AliAnalysisTaskCaloFilter::FillAODv0s()
+{
+  // Copy v0s (use if you know what you do, use quite a lot of memory)
+  
+  if(!fFillv0s) return;
+  
+  // Copy from AODs
+  if(fAODEvent)
+  {
+    TClonesArray* inv0 = fAODEvent  ->GetV0s();
+    TClonesArray* ouv0 =  AODEvent()->GetV0s();
+    
+               //new (ouv0s) TClonesArray(*inv0s);
+    
+    Int_t allv0s = inv0->GetEntriesFast();
+    
+    for (Int_t nv0s = 0; nv0s < allv0s; ++nv0s) 
+    {
+      AliAODv0 *v0 = (AliAODv0*)inv0->At(nv0s);
+      
+      new((*ouv0)[nv0s]) AliAODv0(*v0);
+    } 
     
     return;
   }
   
 }  
 
+//____________________________________________
+void AliAnalysisTaskCaloFilter::FillAODVZERO()
+{
+  // Copy VZERO
+  
+  if(!fFillVZERO) return;
+  
+  AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();
+  
+  if(fESDEvent) *vzeroData = *(fESDEvent->GetVZEROData());
+  else          *vzeroData = *(fAODEvent->GetVZEROData());
+  
+}  
+
 //_______________________________________________
 void AliAnalysisTaskCaloFilter::FillAODVertices()
 {
   // Copy vertices
   
-  AliVEvent* event = InputEvent();
-  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);  
-  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
-  
   // set arrays and pointers
   Double_t pos[3]   ;
   Double_t covVtx[6];
   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
   
   // Copy from AODs
-  if(aodevent)
+  if(fAODEvent)
   {
-    TClonesArray* inVertices = aodevent  ->GetVertices();
+    TClonesArray* inVertices = fAODEvent  ->GetVertices();
     TClonesArray* ouVertices = AODEvent()->GetVertices();
     
-               new (ouVertices) TClonesArray(*inVertices);
+               //new (ouVertices) TClonesArray(*inVertices);
+    
+    //Keep only the first 3 vertices if not requested
+    Int_t allVertices = inVertices->GetEntriesFast();
+    
+    //printf("n Vertices %d\n",allVertices);
+    
+    if(!fFillAllVertices) 
+    {
+      if(allVertices >  3) allVertices = 3;
+    }
+    
+    //printf("Final n Vertices %d\n",allVertices);
+    
+    for (Int_t nVertices = 0; nVertices < allVertices; ++nVertices) 
+    {
+      AliAODVertex *vertex = (AliAODVertex*)inVertices->At(nVertices);
+      
+      new((*ouVertices)[nVertices]) AliAODVertex(*vertex);
+    } 
+    
     return;
   }
   
-  if(!esdevent) return;
+  if(!fESDEvent) return;
   
   // Copy from ESDs
   
@@ -524,169 +880,57 @@ void AliAnalysisTaskCaloFilter::FillAODVertices()
   TClonesArray &vertices = *(AODEvent()->GetVertices());
   
   // Add primary vertex. The primary tracks will be defined
-  // after the loops on the composite objects (V0, cascades, kinks)
-  event->GetPrimaryVertex()->GetXYZ(pos);
-  esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
-  Float_t chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
+  // after the loops on the composite objects (v0, cascades, kinks)
+  fEvent   ->GetPrimaryVertex()->GetXYZ(pos);
+  fESDEvent->GetPrimaryVertex()->GetCovMatrix(covVtx);
+  Float_t chi = fESDEvent->GetPrimaryVertex()->GetChi2toNDF();
   
   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 AliAnalysisTaskCaloFilter::CorrectionsInEMCAL()
+//____________________________________
+void AliAnalysisTaskCaloFilter::Init()
 {
-  //If EMCAL, and requested, correct energy, position ...
-  
-  //Need to do this in a separate loop before filling the ESDs because of the track matching recalculations
+  //Init analysis with configuration macro if available
   
-  if(fCorrect && (fCaloFilter==kEMCAL || fCaloFilter==kBoth) ) 
+  if(gROOT->LoadMacro(fConfigName) >=0)
   {
-    if(!fGeoMatrixSet)
-    {
-      if(fLoadEMCALMatrices)
-      {
-        printf("AliAnalysisTaskCaloFilter::UserExec() - Load user defined EMCAL geometry matrices\n");
-        for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
-          if(fEMCALMatrix[mod]){
-            if(DebugLevel() > 1) 
-              fEMCALMatrix[mod]->Print();
-            fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;  
-          }
-          fGeoMatrixSet=kTRUE;
-        }//SM loop
-      }//Load matrices
-      else if(!gGeoManager)
-      {
-        printf("AliAnalysisTaskCaloFilter::UserExec() - Get geo matrices from data\n");
-        //Still not implemented in AOD, just a workaround to be able to work at least with ESDs        
-        if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) 
-        {
-          if(DebugLevel() > 1) 
-            printf("AliAnalysisTaskCaloFilter Use ideal geometry, values geometry matrix not kept in AODs.\n");
-        }//AOD
-        else 
-        {      
-          if(DebugLevel() > 1) printf("AliAnalysisTaskCaloFilter Load Misaligned matrices. \n");
-          AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
-          if(!esd) 
-          {
-            printf("AliAnalysisTaskCaloFilter::UserExec() - This event does not contain ESDs?");
-            return;
-          }
-          for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
-          {
-            //if(DebugLevel() > 1) 
-            esd->GetEMCALMatrix(mod)->Print();
-            if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
-          } 
-          fGeoMatrixSet=kTRUE;
-        }//ESD
-      }//Load matrices from Data 
-      
-    }//first event
-    
-    //Cluster Loop
-    Int_t nCaloClus = InputEvent()->GetNumberOfCaloClusters();
-    
-    for (Int_t iClust=0; iClust<nCaloClus; ++iClust) 
-    {
-      AliVCluster * cluster = InputEvent()->GetCaloCluster(iClust);
-      
-      if(cluster->IsPHOS()) continue ;
-      
-      Float_t position[]={0,0,0};
-      if(DebugLevel() > 2)
-        printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
-      if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;       
-      
-      if(DebugLevel() > 2)
-      { 
-        printf("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f\n",iClust,cluster->E(),
-               cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel());
-        cluster->GetPosition(position);
-        printf("Filter, before  : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
-      }
-      
-      //Recalculate distance to bad channels, if new list of bad channels provided
-      fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, InputEvent()->GetEMCALCells(), cluster);
-      
-      if(fEMCALRecoUtils->IsRecalibrationOn()) 
-      {
-        fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, InputEvent()->GetEMCALCells());
-        fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
-        fEMCALRecoUtils->RecalculateClusterPID(cluster);
-      }
-      
-      fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
-      
-      if(DebugLevel() > 2)
-      { 
-        printf("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f\n",cluster->GetID(),cluster->E(),
-               cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel());
-        cluster->GetPosition(position);
-        printf("Filter, after   : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
-      }    
-      
-      cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
-      
-    }
+    AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
     
-    //Recalculate track-matching
-    fEMCALRecoUtils->FindMatches(InputEvent(),0,fEMCALGeo);
+    AliAnalysisTaskCaloFilter *filter = (AliAnalysisTaskCaloFilter*)gInterpreter->ProcessLine("ConfigCaloFilter()");
     
-  } // corrections in EMCAL
-}
+    fEMCALGeoName      = filter->fEMCALGeoName; 
+    fLoadEMCALMatrices = filter->fLoadEMCALMatrices;
+    fFillAODFile       = filter->fFillAODFile;
+    fFillTracks        = filter->fFillTracks;
+    fFillHybridTracks  = filter->fFillHybridTracks;
+    fFillv0s           = filter->fFillv0s;
+    fFillVZERO         = filter->fFillVZERO;
+    fFillAllVertices   = filter->fFillAllVertices;
+    fEMCALRecoUtils    = filter->fEMCALRecoUtils; 
+    fConfigName        = filter->fConfigName;
+    fCaloFilter        = filter->fCaloFilter;
+    fEventSelection[0] = filter->fEventSelection[0];
+    fEventSelection[1] = filter->fEventSelection[1];
+    fEventSelection[2] = filter->fEventSelection[2];
+    fAcceptAllMBEvent  = filter->fAcceptAllMBEvent;
+    fMBTriggerMask     = filter->fMBTriggerMask;
+    fCorrect           = filter->fCorrect;
+    fEMCALEnergyCut    = filter->fEMCALEnergyCut;
+    fEMCALNcellsCut    = filter->fEMCALNcellsCut;
+    fPHOSEnergyCut     = filter->fPHOSEnergyCut;
+    fPHOSNcellsCut     = filter->fPHOSNcellsCut;
+    fTrackPtCut        = filter->fTrackPtCut;
+    fVzCut             = filter->fVzCut;
+    fCheckEventVertex  = filter->fCheckEventVertex;
 
-//_______________________________________________________
-Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex()
-{
-  //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
-  //It only works for ESDs
-  
-  AliESDEvent * event = dynamic_cast<AliESDEvent*> (InputEvent());
-  
-  // AODs
-  if(!event) 
-  {
-    // Check that the vertex is not (0,0,0)
-    Double_t v[3];
-    InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
-    
-    if(TMath::Abs(v[2]) < 1e-6 && 
-       TMath::Abs(v[1]) < 1e-6 && 
-       TMath::Abs(v[0]) < 1e-6 ) return kFALSE ;
-    
-    return kTRUE;
-  }
-  
-  // ESDs
-  if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
-  {
-    return kTRUE;
-  }
-  
-  if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) 
-  {
-    // SPD vertex
-    if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) 
-    {
-      //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
-      return kTRUE;
-      
-    }
-    if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) 
-    {
-      //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
-      return kFALSE;
-    }
+    for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
   }
-  return kFALSE;
-  
-}
+} 
 
 //_________________________________________
 void AliAnalysisTaskCaloFilter::PrintInfo()
@@ -694,12 +938,108 @@ void AliAnalysisTaskCaloFilter::PrintInfo()
   //Print settings
   
   printf("AnalysisCaloFilter::PrintInfo() \n");
+  
   printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
   printf("\t Calorimeter Filtering Option     ? %d\n",fCaloFilter);
+  
   //printf("\t Use handmade geo matrices?   EMCAL %d, PHOS %d\n",fLoadEMCALMatrices, fLoadPHOSMatrices);
   printf("\t Use handmade geo matrices?   EMCAL %d, PHOS 0\n",fLoadEMCALMatrices);
-  printf("\t Fill AOD file? %d\n", fFillAODFile);
-  printf("\t Fill Tracks? %d\n"  , fFillTracks);
-  printf("\t Event Selection : EMCAL E > %2.2f, EMCAL nCells >= %d, |vz| < %2.2f\n",fEnergyCut,fNcellsCut,fVzCut);
+  
+  printf("\t Fill: AOD file? %d Tracks? %d; all Vertex? %d; v0s? %d; VZERO ? %d\n", 
+         fFillAODFile,fFillTracks,fFillAllVertices, fFillv0s, fFillVZERO);
+  
+  printf("\t Event Selection based : EMCAL?  %d, PHOS?  %d Tracks? %d - Accept all MB with mask %d? %d\n",
+         fEventSelection[0],fEventSelection[1],fEventSelection[2],fMBTriggerMask, fAcceptAllMBEvent);
+  
+  printf("\t \t EMCAL E > %2.2f, EMCAL nCells >= %d, PHOS E > %2.2f, PHOS nCells >= %d, Track pT > %2.2f, |vz| < %2.2f\n",
+         fEMCALEnergyCut,fEMCALNcellsCut,fPHOSEnergyCut,fPHOSNcellsCut, fTrackPtCut,fVzCut);
+}
+
+//_______________________________________________________
+void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
+{
+  // Init geometry 
+       
+  fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;  
+  
+  if(fFillMCParticles)
+  {
+    TClonesArray * aodMCParticles = new TClonesArray("AliAODMCParticle",500);
+               aodMCParticles->SetName("mcparticles");
+               ((AliAODHandler*)AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler())->AddBranch("TClonesArray", &aodMCParticles);
+  }
+  
+}  
+
+//____________________________________________________________
+void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
+{
+  // Execute analysis for current event
+  // Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
+  
+  AliDebug(1,Form("Analysing event # %d", (Int_t)Entry()));
+  
+  fEvent    = InputEvent();
+  fAODEvent = dynamic_cast<AliAODEvent*> (fEvent);  
+  fESDEvent = dynamic_cast<AliESDEvent*> (fEvent);
+  
+  if(!fEvent) 
+  {
+    AliInfo("This event does not contain Input?");
+    return;
+  }
+  
+  // printf("Start processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
+  
+  // Select the event
+  
+  if(!AcceptEvent()) return ;
+  
+  //Magic line to write events to file
+  
+  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
+  
+  // Reset output AOD
+  
+  Int_t nVertices = 0;
+  if(fFillv0s) nVertices = fEvent->GetNumberOfV0s();
+  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
+  Int_t nTracks   = fEvent->GetNumberOfTracks();
+  
+  AODEvent()->ResetStd(nTracks, nVertices, 0, 0, 0, nCaloClus, 0, 0);
+  
+  // Copy
+  
+  FillAODHeader();
+  
+  // 
+  FillAODv0s();
+  
+  //
+  FillAODVertices(); // Do it before the track filtering to have the reference to the vertex
+  
+  // 
+  FillAODVZERO();
+  
+  //
+  FillAODTracks();
+  
+  //
+  CorrectionsInEMCAL();
+  
+  //
+  FillAODCaloClusters();
+  
+  //
+  FillAODCaloCells();
+  
+  //
+  FillAODCaloTrigger();
+  
+  // 
+  FillAODMCParticles();
+  
+  //printf("Filtered event, end processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
+  
 }