added copy of VZERO, copy only hybrid tracks, copy only the main vertices, select...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 May 2012 16:29:22 +0000 (16:29 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 May 2012 16:29:22 +0000 (16:29 +0000)
PWGGA/CaloTasks/AliAnalysisTaskCaloFilter.cxx
PWGGA/CaloTasks/AliAnalysisTaskCaloFilter.h
PWGGA/CaloTasks/macros/AddTaskCaloFilter.C

index d844b55..20a9770 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)
 //////////////////////////////////////////////////////////
@@ -50,17 +51,29 @@ ClassImp(AliAnalysisTaskCaloFilter)
 ////////////////////////////////////////////////////////////////////////
 
 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
-AliAnalysisTaskSE("CaloFilterTask"), //fCuts(0x0),
-fCaloFilter(0), fCorrect(kFALSE), 
-fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"), 
+AliAnalysisTaskSE("CaloFilterTask"), 
+fCaloFilter(0),           fEventSelection(), 
+fAcceptAllMBEvent(kFALSE),  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), 
+fFillTracks(kFALSE),      fFillHybridTracks(kFALSE),
+fFillAllVertices(kFALSE), fFillv0s(kFALSE),  
+fFillVZERO(kFALSE),
+fEMCALEnergyCut(0.),      fEMCALNcellsCut (0),
+fPHOSEnergyCut(0.),       fPHOSNcellsCut (0), 
+fTrackPtCut(-1),
+fVzCut(100.),             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 ;
   
@@ -68,17 +81,29 @@ 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),  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),
+fFillTracks(kFALSE),      fFillHybridTracks(kFALSE),
+fFillAllVertices(kFALSE), fFillv0s(kFALSE),
+fFillVZERO(kFALSE),
+fEMCALEnergyCut(0.),      fEMCALNcellsCut(0), 
+fPHOSEnergyCut(0.),       fPHOSNcellsCut(0), 
+fTrackPtCut(-1),
+fVzCut(100.),             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 ;
   
@@ -94,91 +119,129 @@ 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;
-    
-    for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
-  }
-} 
+  Bool_t eventSel = kFALSE;
+  
+  Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+  
+  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
+  
+  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))
+    { 
+      
+      if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventEMCAL() - Accept :  E %2.2f > %2.2f, nCells %d > %d \n", 
+                             clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut);
+      
+      return kTRUE;
+    }
+    
+  }// loop
+  
+  if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventEMCAL() - Reject \n");
+  
+  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 ))
+    { 
+      
+      if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventPHOS() - Accept :  E %2.2f > %2.2f, nCells %d > %d \n", 
+                             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();
+  if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventPHOS() - Reject \n");
   
-  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) 
+    {
+      if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventTrack() - Accept :  pT %2.2f > %2.2f \n", 
+                             momentum.Pt(), fTrackPtCut);
+
+      return kTRUE;
+    }
+    
+  } 
   
-  //
-  FillAODCaloClusters();
+  if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventTrack() - Reject \n");
   
-  //
-  FillAODCaloCells();
+  return kFALSE;
   
-}
+}  
 
 //___________________________________________________
 Bool_t AliAnalysisTaskCaloFilter::AcceptEventVertex()
@@ -188,53 +251,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) 
+  {
+    if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventVertex() - Vz Reject \n");
+    return kFALSE ;
+  }
   
   return CheckForPrimaryVertex();
 }  
 
-//__________________________________________________
-Bool_t AliAnalysisTaskCaloFilter::AcceptEventEMCAL()
+//_______________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex()
 {
-  // Accept event given there is a cluster with enough energy
-  
-  if(fCaloFilter!=kEMCAL) return kTRUE; // accept all
-  
-  Int_t           nCluster = InputEvent() -> GetNumberOfCaloClusters();
-  AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
-  Int_t           bc       = InputEvent() -> GetBunchCrossNumber();
+  //Check if the vertex was well reconstructed, copy from v0Reader of conversion group
+  //It only works for ESDs
   
-  for(Int_t icalo = 0; icalo < nCluster; icalo++)
+  // AODs
+  if(!fESDEvent) 
   {
-    AliESDCaloCluster *clus = (AliESDCaloCluster*) (InputEvent()->GetCaloCluster(icalo));
+    // Check that the vertex is not (0,0,0)
+    Double_t v[3];
+    InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
     
-    if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fNcellsCut ) && ( clus->E() > fEnergyCut ) &&
-       fEMCALRecoUtils->IsGoodCluster(clus,fEMCALGeo,caloCell,bc))
-    { 
-      
-      // printf("Accept event %d, E %2.2f > %2.2f, nCells %d > %d \n", 
-      //        (Int_t) Entry(),clus->E(), fEnergyCut, clus->GetNCells(), fNcellsCut);
+    if(TMath::Abs(v[2]) < 1e-6 && 
+       TMath::Abs(v[1]) < 1e-6 && 
+       TMath::Abs(v[0]) < 1e-6 ) 
+    {
+      if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject \n");
       
-      return kTRUE;
+      return kFALSE ;
     }
     
-  }// loop
+    return kTRUE;
+  }
+  
+  // ESDs
+  if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() > 0)
+  {
+    return kTRUE;
+  }
+  
+  if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() < 1) 
+  {
+    // SPD vertex
+    if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() > 0) 
+    {
+      //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+      return kTRUE;
+      
+    }
+    if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() < 1) 
+    {
+      //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+      if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject \n");
+      return kFALSE;
+    }
+  }
+  
+  if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject \n");
   
   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)
+      {
+        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));
+      
+    }
+    
+    //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());
@@ -268,9 +461,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());
@@ -293,19 +486,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.
     
@@ -392,30 +583,45 @@ void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
   
 }
 
+//__________________________________________________
+void AliAnalysisTaskCaloFilter::FillAODCaloTrigger()
+{
+  // AOD CaloTrigger copy
+  
+  AliAODCaloTrigger* triggerEM = AODEvent()->GetCaloTrigger("EMCAL");
+  AliAODCaloTrigger* triggerPH = AODEvent()->GetCaloTrigger("PHOS");
+  
+  // Copy from AODs
+  if(fAODEvent)
+  {
+    if(fCaloFilter==kBoth || fCaloFilter==kPHOS)  *triggerPH = *(fAODEvent->GetCaloTrigger("PHOS"));
+    
+    if(fCaloFilter==kBoth || fCaloFilter==kEMCAL) *triggerEM = *(fAODEvent->GetCaloTrigger("EMCAL"));
+    
+    return;
+  }
+  
+}  
 
 //______________________________________________
 void AliAnalysisTaskCaloFilter::FillAODHeader()
 {
   // AOD header copy
   
-  AliVEvent* event = InputEvent();
-  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);  
-  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
-  
   AliAODHeader* header = AODEvent()->GetHeader();
   
   // Copy from AODs
-  if(aodevent)
+  if(fAODEvent)
   {
-    *header = *(aodevent->GetHeader());
+    *header = *(fAODEvent->GetHeader());
     return;
   }
   
-  if(!esdevent) return;
+  if(!fESDEvent) return;
   
   // Copy from ESDs
   
-  header->SetRunNumber(event->GetRunNumber());
+  header->SetRunNumber(fEvent->GetRunNumber());
   
   TTree* tree = fInputHandler->GetTree();
   if (tree) 
@@ -424,15 +630,15 @@ void AliAnalysisTaskCaloFilter::FillAODHeader()
     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
   {
@@ -441,27 +647,27 @@ 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]={fEvent->GetDiamondX(),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());
   
 }
 
@@ -472,52 +678,168 @@ 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);
-        
-    for (Int_t nTrack = 0; nTrack < aodevent->GetNumberOfTracks(); ++nTrack) 
+    
+    //printf("N tracks %d\n",fAODEvent->GetNumberOfTracks());
+    Int_t nCopyTrack = 0;
+    for (Int_t nTrack = 0; nTrack < fAODEvent->GetNumberOfTracks(); ++nTrack) 
     {
-      AliAODTrack *track = aodevent->GetTrack(nTrack);
+      AliAODTrack *track = fAODEvent->GetTrack(nTrack);
+      
+      // 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!!!
 
-      new((*ouTracks)[nTrack]) AliAODTrack(*track);
+      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(),
+                                                            track->Chi2perNDF());
+      
+      
+      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->SetTrackPhiEtaOnEMCal( track->GetTrackPhiOnEMCal(),  track->GetTrackPhiOnEMCal() );
+          
     } 
     
+    //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);
     
-    for (Int_t nVertices = 0; nVertices < inVertices->GetEntriesFast(); ++nVertices) 
+    //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);
       
@@ -527,7 +849,7 @@ void AliAnalysisTaskCaloFilter::FillAODVertices()
     return;
   }
   
-  if(!esdevent) return;
+  if(!fESDEvent) return;
   
   // Copy from ESDs
   
@@ -536,169 +858,55 @@ 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));
-      
-    }
-    
-    //Recalculate track-matching
-    fEMCALRecoUtils->FindMatches(InputEvent(),0,fEMCALGeo);
+    printf("Configure analysis with %s\n",fConfigName.Data());
     
-  } // corrections in EMCAL
-}
-
-//_______________________________________________________
-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) ;
+    AliAnalysisTaskCaloFilter *filter = (AliAnalysisTaskCaloFilter*)gInterpreter->ProcessLine("ConfigCaloFilter()");
     
-    if(TMath::Abs(v[2]) < 1e-6 && 
-       TMath::Abs(v[1]) < 1e-6 && 
-       TMath::Abs(v[0]) < 1e-6 ) return kFALSE ;
+    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;
+    fCorrect           = filter->fCorrect;
+    fEMCALEnergyCut    = filter->fEMCALEnergyCut;
+    fEMCALNcellsCut    = filter->fEMCALNcellsCut;
+    fPHOSEnergyCut     = filter->fPHOSEnergyCut;
+    fPHOSNcellsCut     = filter->fPHOSNcellsCut;
+    fTrackPtCut        = filter->fTrackPtCut;
+    fVzCut             = filter->fVzCut;
     
-    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()
@@ -706,12 +914,99 @@ 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? %d\n",
+         fEventSelection[0],fEventSelection[1],fEventSelection[2],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) ;  
+  
+}  
+
+//____________________________________________________________
+void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
+{
+  // 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());
+  
+  fEvent    = InputEvent();
+  fAODEvent = dynamic_cast<AliAODEvent*> (fEvent);  
+  fESDEvent = dynamic_cast<AliESDEvent*> (fEvent);
+  
+  if(!fEvent) 
+  {
+    printf("AliAnalysisTaskCaloFilter::UserExec - 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();
+  
+  //printf("Filtered event, end processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
+  
 }
 
index 8f999f8..f75055e 100644 (file)
@@ -8,8 +8,8 @@
 // Filter the ESDCaloClusters and ESDCaloCells of EMCAL,
 // PHOS or both, creating the corresponing AODCaloClusters
 // and AODCaloCells.
-// Keep also the AODHeader information and the vertex.
-// Keep tracks if requested
+// Fill also the AODHeader information and the vertex.
+// Fill tracks if requested
 // Copy of AliAnalysisTaskESDfilter.
 // Author: Gustavo Conesa Balbastre (INFN - Frascati)
 //////////////////////////////////////////////////////////
 #include "AliAnalysisTaskSE.h"
 class AliEMCALRecoUtils;
 class AliEMCALGeometry;
+class AliESDEvent;
+class AliAODEvent;
 
 class AliAnalysisTaskCaloFilter : public AliAnalysisTaskSE
 {
  public:
+  
   AliAnalysisTaskCaloFilter();
   AliAnalysisTaskCaloFilter(const char* name);
   virtual ~AliAnalysisTaskCaloFilter() ;
@@ -32,7 +35,15 @@ class AliAnalysisTaskCaloFilter : public AliAnalysisTaskSE
   virtual void   LocalInit() { Init() ; }
   virtual void   UserExec(Option_t *option);
         
+  // Task own methods
+  
+  Bool_t  AcceptEvent() ;
+  
   Bool_t  AcceptEventEMCAL();
+
+  Bool_t  AcceptEventPHOS();
+  
+  Bool_t  AcceptEventTrack();
   
   Bool_t  AcceptEventVertex();
   
@@ -46,10 +57,16 @@ class AliAnalysisTaskCaloFilter : public AliAnalysisTaskSE
   
   void    FillAODCaloClusters();
   
+  void    FillAODCaloTrigger();
+
   void    FillAODTracks();
+
+  void    FillAODv0s();
   
   void    FillAODVertices();
   
+  void    FillAODVZERO();
+  
   void    PrintInfo();
   
   // * Task settings *
@@ -67,35 +84,67 @@ class AliAnalysisTaskCaloFilter : public AliAnalysisTaskSE
   //void    SwitchOffLoadOwnPHOSGeometryMatrices()              { fLoadPHOSMatrices = kFALSE ; }
   //void    SetPHOSGeometryMatrixInSM(TGeoHMatrix* m, Int_t i)  { fPHOSMatrix[i]    = m      ; }
     
-  void    SwitchOnFillAODFile()                   { fFillAODFile = kTRUE    ; }
-  void    SwitchOffFillAODFile()                  { fFillAODFile = kFALSE   ; }
+  void    SwitchOnFillAODFile()                   { fFillAODFile = kTRUE        ; }
+  void    SwitchOffFillAODFile()                  { fFillAODFile = kFALSE       ; }
 
-  void    SwitchOnFillTracks()                    { fFillTracks  = kTRUE    ; }
-  void    SwitchOffFillTracks()                   { fFillTracks  = kFALSE   ; }
+  void    SwitchOnFillTracks()                    { fFillTracks  = kTRUE        ; }
+  void    SwitchOffFillTracks()                   { fFillTracks  = kFALSE       ; }
+  
+  void    SwitchOnFillHybridTracks()              { fFillTracks  = kTRUE        ;
+                                                    fFillHybridTracks  = kTRUE  ; }
+  void    SwitchOffFillHybridTracks()             { fFillHybridTracks  = kFALSE ; }
+  
+  void    SwitchOnFillv0s()                       { fFillv0s     = kTRUE        ; }
+  void    SwitchOffFillv0s()                      { fFillv0s     = kFALSE       ; }
+  
+  void    SwitchOnFillVZERO()                     { fFillVZERO   = kTRUE        ; }
+  void    SwitchOffFillVZERO()                    { fFillVZERO   = kFALSE       ; }
+  
+  void    SwitchOnFillAllVertices()               { fFillAllVertices = kTRUE    ; }
+  void    SwitchOffFillAllVertices()              { fFillAllVertices = kFALSE   ; }
   
   enum    caloFilter {kBoth = 0, kEMCAL = 1, kPHOS=2};
-  void    SetCaloFilter(Int_t calo)               { fCaloFilter = calo      ; }
-  TString GetCaloFilter()                  const  { return fCaloFilter      ; }  
   
-  void    SetEMCALRecoUtils(AliEMCALRecoUtils* ru){ fEMCALRecoUtils = ru    ; }
-  AliEMCALRecoUtils* GetEMCALRecoUtils()   const  { return fEMCALRecoUtils  ; }
+  void    SetCaloFilter(Int_t calo)               { fCaloFilter = calo          ; }
+  TString GetCaloFilter()                  const  { return fCaloFilter          ; }  
+  
+  void    SetEventSelection(Bool_t emcal, Bool_t phos, Bool_t track) { 
+    // Detector involved in event selection
+    fEventSelection[0] = emcal ; fEventSelection[1] = phos ; fEventSelection[2] = track ; }
+  
+  void    SwitchOnAcceptAllMBEvent()              { fAcceptAllMBEvent = kTRUE   ; }
+  void    SwitchOffAcceptAllMBEvent()             { fAcceptAllMBEvent = kFALSE  ; }
+
+  void    SetEMCALRecoUtils(AliEMCALRecoUtils* ru){ fEMCALRecoUtils = ru        ; }
+  AliEMCALRecoUtils* GetEMCALRecoUtils()   const  { return fEMCALRecoUtils      ; }
 
-  void    SwitchOnClusterCorrection()             { fCorrect = kTRUE        ; }
-  void    SwitchOffClusterCorrection()            { fCorrect = kFALSE       ; }
+  void    SwitchOnClusterCorrection()             { fCorrect = kTRUE            ; }
+  void    SwitchOffClusterCorrection()            { fCorrect = kFALSE           ; }
+  
+  void    SetConfigFileName(TString name)         { fConfigName = name          ; }
   
-  void    SetConfigFileName(TString name)         { fConfigName = name      ; }
+  void    SetEMCALEnergyCut(Float_t cut)          { fEMCALEnergyCut = cut       ; }
+  Float_t GetEMCALEnergyCut()               const { return fEMCALEnergyCut      ; }
+  void    SetEMCALNcellsCut(Int_t cut)            { fEMCALNcellsCut = cut       ; }
+  Int_t   GetEMCALNcellsCut()               const { return fEMCALNcellsCut      ; }
   
-  void    SetEnergyCut(Float_t cut)               { fEnergyCut = cut        ; }
-  Float_t GetEnergyCut()                    const { return fEnergyCut       ; }
-  void    SetNcellsCut(Int_t cut)                 { fNcellsCut = cut        ; }
-  Int_t   GetNcellsCut()                    const { return fNcellsCut       ; }
-  void    SetVzCut(Float_t cut)                   { fVzCut = cut            ; }
-  Float_t GetVzCut()                        const { return fVzCut           ; }
+  void    SetPHOSEnergyCut(Float_t cut)           { fPHOSEnergyCut = cut        ; }
+  Float_t GetPHOSEnergyCut()                const { return fPHOSEnergyCut       ; }
+  void    SetPHOSNcellsCut(Int_t cut)             { fPHOSNcellsCut = cut        ; }
+  Int_t   GetPHOSNcellsCut()                const { return fPHOSNcellsCut       ; }
+  
+  void    SetTrackPtCut(Float_t cut)              { fTrackPtCut = cut           ; }
+  Float_t GetTrackPtCut()                   const { return fTrackPtCut          ; }
+  
+  void    SetVzCut(Float_t cut)                   { fVzCut = cut                ; }
+  Float_t GetVzCut()                        const { return fVzCut               ; }
   
   
 private:
     
   Int_t               fCaloFilter;        // Calorimeter to filter
+  Bool_t              fEventSelection[3]; // Define which detector is used to select the event
+  Bool_t              fAcceptAllMBEvent;    // Do not select the MB events with same cuts as other triggers
   Int_t               fCorrect;           // Recalibrate or recalculate different cluster parameters
   
   //EMCAL specific
@@ -111,17 +160,34 @@ private:
   Bool_t              fGeoMatrixSet;      // Set geometry matrices only once, for the first event.   
   
   TString             fConfigName;        // Name of analysis configuration file
+  
   Bool_t              fFillAODFile;       // Fill the output AOD file with clusters 
   Bool_t              fFillTracks;        // Fill tracks
+  Bool_t              fFillHybridTracks;  // Fill hybrid tracks
+
+  Bool_t              fFillAllVertices;   // Fill all vertices
+  Bool_t              fFillv0s;           // Fill v0s
+  Bool_t              fFillVZERO;         // Fill VZERO
+
+  Float_t             fEMCALEnergyCut;    //  At least an EMCAL cluster with this energy in the event
+  Int_t               fEMCALNcellsCut;    //  At least an EMCAL cluster with fNCellsCut cells over fEnergyCut
+
+  Float_t             fPHOSEnergyCut;     //  At least a PHOS cluster with this energy in the event
+  Int_t               fPHOSNcellsCut;     //  At least a PHOS cluster with fNCellsCut cells over fEnergyCut
+  
+  Float_t             fTrackPtCut;        //  At least a track with this pT in the event
   
-  Float_t             fEnergyCut;         //  At least a cluster with this energy in the event
-  Int_t               fNcellsCut;         //  At least a cluster with fNCellsCut cells over fEnergyCut
   Float_t             fVzCut;             //  At least events with vertex within cut
   
+  AliVEvent*          fEvent;             //! event pointer
+  AliESDEvent*        fESDEvent;          //! ESD event pointer
+  AliAODEvent*        fAODEvent;          //! AOD event pointer
+
+  
   AliAnalysisTaskCaloFilter(           const AliAnalysisTaskCaloFilter&);
   AliAnalysisTaskCaloFilter& operator=(const AliAnalysisTaskCaloFilter&);
   
-  ClassDef(AliAnalysisTaskCaloFilter, 8); // Analysis task for standard ESD filtering
+  ClassDef(AliAnalysisTaskCaloFilter, 9); // Analysis task for standard ESD filtering
   
 };
 
index c575886..a1c122f 100644 (file)
@@ -1,8 +1,8 @@
 AliAnalysisTaskCaloFilter * AddTaskCaloFilter(const Bool_t  bias      = kTRUE, 
-                                              const Float_t minE      = 6, 
-                                              const Float_t minN      = 2, 
+                                              const Float_t minE      = 5, 
+                                              const Float_t minN      = 3, 
                                               const Float_t vz        = 10.,
-                                              const Int_t   opt       = AliAnalysisTaskCaloFilter::kEMCAL, //kPHOS, kEMCAL or kBoth
+                                              const Int_t   opt       = AliAnalysisTaskCaloFilter::kBoth, //kPHOS, kEMCAL or kBoth
                                               const Bool_t  correct   = kFALSE,
                                               const Bool_t  fillTrack = kTRUE,
                                               const Bool_t  fillAOD   = kTRUE)
@@ -27,25 +27,44 @@ AliAnalysisTaskCaloFilter * AddTaskCaloFilter(const Bool_t  bias      = kTRUE,
   
   AliAnalysisTaskCaloFilter * filter = new AliAnalysisTaskCaloFilter("CaloFilter");
   
+  //filter->SetDebugLevel(2);
+  
   filter->SetCaloFilter(opt); //kPHOS, kEMCAL or kBoth
   
   filter->SetVzCut(vz);
   
-  if(bias) // select events with significant signal in EMCAL
+  if(bias) // select events with significant signal in EMCAL or TPC or PHOS
   {
-    filter->SetEnergyCut(minE);
-    filter->SetNcellsCut(minN);
+    filter->SetEventSelection(1,1,1); // Select events depending on EMCAL, PHOS and tracks criteria
+    filter->SwitchOnAcceptAllMBEvent();
+    
+    filter->SetEMCALEnergyCut(minE);
+    filter->SetEMCALNcellsCut(minN);
+    filter->SetPHOSEnergyCut(minE);
+    filter->SetPHOSNcellsCut(minN);
     
+    filter->SetTrackPtCut(minE);
+
     filter->SelectCollisionCandidates(AliVEvent::kAny) ;
     
     printf("--- Select events with 1 cluster with at least %2.2f GeV and N = %d ---\n",minE,minN);
   }
   else // Do not bias the signal in EMCAL, select MB events 
   {
-    filter->SetEnergyCut(0);
-    filter->SetNcellsCut(0);    
-
-    filter->SelectCollisionCandidates(AliVEvent::kAnyINT | AliVEvent::kCentral | AliVEvent::kSemiCentral ) ;
+    
+    filter->SetEventSelection(0,0,0);
+    filter->SwitchOnAcceptAllMBEvent();
+    
+    filter->SetEMCALEnergyCut(-1);
+    filter->SetEMCALNcellsCut(0);  
+    
+    filter->SetPHOSEnergyCut(-1);
+    filter->SetPHOSNcellsCut(0); 
+    
+    filter->SetTrackPtCut(-1);
+    
+    filter->SelectCollisionCandidates(AliVEvent::kMB);// | AliVEvent::kCentral | AliVEvent::kSemiCentral ) ;
     
     printf("--- Select Min Bias events ---\n");
   }
@@ -53,8 +72,17 @@ AliAnalysisTaskCaloFilter * AddTaskCaloFilter(const Bool_t  bias      = kTRUE,
   if(correct)   filter->SwitchOnClusterCorrection();
   else          filter->SwitchOffClusterCorrection();  
   
-  if(fillTrack) filter->SwitchOnFillTracks();
-  else          filter->SwitchOffFillTracks();
+  AliEMCALRecoUtils * reco = filter->GetEMCALRecoUtils();
+  reco->SwitchOnRejectExoticCluster() ;
+  reco->SetExoticCellFractionCut(0.97);
+  reco->SetExoticCellMinAmplitudeCut(2.);
+
+  if(fillTrack) { filter->SwitchOnFillTracks()  ; filter->SwitchOnFillHybridTracks()  ; }
+  else          { filter->SwitchOffFillTracks() ; filter->SwitchOffFillHybridTracks() ; }
+  
+  filter->SwitchOffFillv0s() ; // Put ON if you know what you do.
+  
+  filter->SwitchOnFillVZERO(); // Be able to recalculate centrality and event plane afterwards even it is stored in header
   
   if(fillAOD)   filter->SwitchOnFillAODFile();
   else          filter->SwitchOffFillAODFile();