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 d844b5512194eb7866c05b8baa64953bfb528ac0..20a9770bc447788dce57c2ffce14c2d17bf5a6aa 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,111 +119,35 @@ AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
   
 }
 
-//-------------------------------------------------------------------
-void AliAnalysisTaskCaloFilter::Init()
-{
-  //Init analysis with configuration macro if available
-  
-  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] ;
-  }
-} 
-
-//__________________________________________________
-void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
-{
-  // Init geometry 
-       
-  fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;  
-  
-}  
-
-//____________________________________________________________
-void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
+//_____________________________________________
+Bool_t AliAnalysisTaskCaloFilter::AcceptEvent()
 {
-  // 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());
-  
-  AliVEvent* event = InputEvent();
-  if(!event) 
-  {
-    printf("AliAnalysisTaskCaloFilter::UserExec - This event does not contain Input?");
-    return;
-  }
+  // Define conditions to accept the event to be filtered
   
-  // Select the event
-  
-  if(!AcceptEventVertex()) return;
-  if(!AcceptEventEMCAL())  return;
-  
-  //Magic line to write events to file
-  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
-  
-  Int_t nVertices = event->GetNumberOfV0s()+event->GetNumberOfCascades();;
-  Int_t nCaloClus = event->GetNumberOfCaloClusters();
-  Int_t nTracks   = event->GetNumberOfTracks();
-  
-  AODEvent()->ResetStd(nTracks, nVertices, 0, 0, 0, nCaloClus, 0, 0);
+  Bool_t eventSel = kFALSE;
   
-  //
-  FillAODHeader();
+  Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
   
-  //
-  FillAODVertices();
+  if     ( isMB && fAcceptAllMBEvent )                eventSel = kTRUE; // accept any MB event
   
-  //
-  FillAODTracks();
+  else if( fEventSelection[0] && AcceptEventEMCAL() ) eventSel = kTRUE; // accept event depending on EMCAL activity
   
-  //
-  CorrectionsInEMCAL();
+  else if( fEventSelection[1] && AcceptEventPHOS () ) eventSel = kTRUE; // accept event depending on PHOS  activity
   
-  //
-  FillAODCaloClusters();
+  else if( fEventSelection[2] && AcceptEventTrack() ) eventSel = kTRUE; // accept event depending on Track activity
   
-  //
-  FillAODCaloCells();
+  return eventSel ;
   
 }
 
-//___________________________________________________
-Bool_t AliAnalysisTaskCaloFilter::AcceptEventVertex()
-{
-  // Accept event with good vertex
-  
-  Double_t v[3];
-  InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
-  
-  if(TMath::Abs(v[2]) > fVzCut) return kFALSE ;
-  
-  return CheckForPrimaryVertex();
-}  
-
 //__________________________________________________
 Bool_t AliAnalysisTaskCaloFilter::AcceptEventEMCAL()
 {
-  // Accept event given there is a cluster with enough energy
+  // Accept event given there is a EMCAL cluster with enough energy, and not noisy, exotic
+  
+  if(fCaloFilter==kPHOS)   return kTRUE; // accept 
   
-  if(fCaloFilter!=kEMCAL) return kTRUE; // accept all
+  if(fEMCALEnergyCut <= 0) return kTRUE; // accept
   
   Int_t           nCluster = InputEvent() -> GetNumberOfCaloClusters();
   AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
@@ -206,345 +155,162 @@ Bool_t AliAnalysisTaskCaloFilter::AcceptEventEMCAL()
   
   for(Int_t icalo = 0; icalo < nCluster; icalo++)
   {
-    AliESDCaloCluster *clus = (AliESDCaloCluster*) (InputEvent()->GetCaloCluster(icalo));
+    AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
     
-    if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fNcellsCut ) && ( clus->E() > fEnergyCut ) &&
+    if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
        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 (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::FillAODCaloCells()
-{  
-  // Fill EMCAL/PHOS cell info
+//_________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::AcceptEventPHOS()
+{
+  // Accept event given there is a PHOS cluster with enough energy and not noisy/exotic
   
-  AliVEvent * event = InputEvent();
+  if(fCaloFilter==kEMCAL) return kTRUE; // accept 
   
-  // EMCAL
-  if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && event->GetEMCALCells()) 
-  { // protection against missing ESD information
-    AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
-    Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
+  if(fPHOSEnergyCut <= 0) return kTRUE; // accept
+  
+  Int_t nCluster = InputEvent() -> GetNumberOfCaloClusters();
+  
+  for(Int_t icalo = 0; icalo < nCluster; icalo++)
+  {
+    AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
     
-    AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
-    aodEMcells.CreateContainer(nEMcell);
-    aodEMcells.SetType(AliVCaloCells::kEMCALCell);
-    Double_t calibFactor = 1.;   
-    for (Int_t iCell = 0; iCell < nEMcell; iCell++) 
+    if( ( clus->IsPHOS() ) && ( clus->GetNCells() > fPHOSNcellsCut ) && ( clus->E() > fPHOSEnergyCut ))
     { 
-      Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
-      fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); 
-      fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);      
       
-      if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn())
-      { 
-        calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
-      }
+      if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventPHOS() - Accept :  E %2.2f > %2.2f, nCells %d > %d \n", 
+                             clus->E(), fPHOSEnergyCut, clus->GetNCells(), fPHOSNcellsCut);
       
-      if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
-      { //Channel is not declared as bad
-        aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
-                           eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
-        //printf("GOOD channel\n");
-      }
-      else 
-      {
-        aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
-        //printf("BAD channel\n");
-      }
-    }
-    aodEMcells.Sort();
-  }
-  
-  // PHOS 
-  if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && event->GetPHOSCells()) 
-  { // protection against missing ESD information
-    AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
-    Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
-    
-    AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
-    aodPHcells.CreateContainer(nPHcell);
-    aodPHcells.SetType(AliVCaloCells::kPHOSCell);
-    
-    for (Int_t iCell = 0; iCell < nPHcell; iCell++) 
-    {      
-      aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell),
-                         eventPHcells.GetTime(iCell),eventPHcells.GetMCLabel(iCell),eventPHcells.GetEFraction(iCell));
+      return kTRUE;
     }
     
-    aodPHcells.Sort();
-  }
-}
-
+  }// loop
+  
+  if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventPHOS() - Reject \n");
+  
+  return kFALSE;
+  
+}  
 
-//___________________________________________________
-void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
+//__________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::AcceptEventTrack()
 {
-  // Fill the AOD with caloclusters
-  
-  AliVEvent * event = InputEvent();
+  // Accept event if there is a track avobe a certain pT
   
-  // Access to the AOD container of clusters
+  if(fTrackPtCut <= 0) return kTRUE; // accept
   
-  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
-  Int_t jClusters=0;
-  Float_t  posF[3]  ;
+  Double_t pTrack[3] = {0,0,0};
   
-  Int_t nCaloClus = event->GetNumberOfCaloClusters();
-  for (Int_t iClust=0; iClust<nCaloClus; ++iClust) 
+  for (Int_t nTrack = 0; nTrack < fEvent->GetNumberOfTracks(); ++nTrack) 
   {
+    AliVTrack *track = (AliVTrack*) fEvent->GetTrack(nTrack);
     
-    AliVCluster * cluster = event->GetCaloCluster(iClust);
-    
-    //Check which calorimeter information we want to keep.
-    
-    if(fCaloFilter!=kBoth)
-    {
-      if     (fCaloFilter==kPHOS  && cluster->IsEMCAL()) continue ;
-      else if(fCaloFilter==kEMCAL && cluster->IsPHOS())  continue ;
-    }  
-    
-    // Get original residuals, in case of previous recalculation, reset them 
-    Float_t dR = cluster->GetTrackDx();
-    Float_t dZ = cluster->GetTrackDz();
+    // Select only hybrid tracks?
+    if(fAODEvent && fFillHybridTracks && !((AliAODTrack*)track)->IsHybridGlobalConstrainedGlobal()) continue;
     
-    if(DebugLevel() > 2)
-      printf("Original residuals : dZ %f, dR %f\n ",dZ, dR);
+    track->GetPxPyPz(pTrack) ;
     
-    //--------------------------------------------------------------
-    //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());
-      
-      if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;       
-      
-      if(fEMCALRecoUtils->IsExoticCluster(cluster, InputEvent()->GetEMCALCells(),InputEvent()->GetBunchCrossNumber())) continue;       
-      
-      fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
-      cluster->SetTrackDistance(dR,dZ);
-    }
+    TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
     
-    if(DebugLevel() > 2)
+    if(momentum.Pt() > fTrackPtCut) 
     {
-      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);
+      if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventTrack() - Accept :  pT %2.2f > %2.2f \n", 
+                             momentum.Pt(), fTrackPtCut);
+
+      return kTRUE;
     }
     
-    //--------------------------------------------------------------
-    
-    //Now fill AODs
-    
-    Int_t id       = cluster->GetID();
-    Float_t energy = cluster->E();
-    cluster->GetPosition(posF);
-    
-    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
-    AliAODCaloCluster(id,
-                      cluster->GetNLabels(),
-                      cluster->GetLabels(),
-                      energy,
-                      posF,
-                      NULL,
-                      cluster->GetType());
-    
-    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
-                                cluster->GetDispersion(),
-                                cluster->GetM20(), cluster->GetM02(),
-                                -1,  
-                                cluster->GetNExMax(),cluster->GetTOF()) ;
-    
-    caloCluster->SetPIDFromESD(cluster->GetPID());
-    caloCluster->SetNCells(cluster->GetNCells());
-    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
-    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]);
-    }    
-    
-    //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) 
-    { //Default value in PHOS 999, in EMCAL 1024, why?
-      caloCluster->AddTrackMatched(new AliAODTrack); 
-    }
-    // TO DO, in case Tracks available, think how to put the matched track in AOD
-  }
+  } 
   
-  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots  
+  if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventTrack() - Reject \n");
   
-}
-
+  return kFALSE;
+  
+}  
 
-//______________________________________________
-void AliAnalysisTaskCaloFilter::FillAODHeader()
+//___________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::AcceptEventVertex()
 {
-  // AOD header copy
-  
-  AliVEvent* event = InputEvent();
-  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);  
-  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
+  // Accept event with good vertex
   
-  AliAODHeader* header = AODEvent()->GetHeader();
+  Double_t v[3];
+  InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
   
-  // Copy from AODs
-  if(aodevent)
+  if(TMath::Abs(v[2]) > fVzCut) 
   {
-    *header = *(aodevent->GetHeader());
-    return;
+    if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventVertex() - Vz Reject \n");
+    return kFALSE ;
   }
   
-  if(!esdevent) return;
+  return CheckForPrimaryVertex();
+}  
+
+//_______________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex()
+{
+  //Check if the vertex was well reconstructed, copy from v0Reader of conversion group
+  //It only works for ESDs
   
-  // Copy from ESDs
-  
-  header->SetRunNumber(event->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());
-  
-  //Centrality
-  if(event->GetCentrality())
-  {
-    header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
-  }
-  else
-  {
-    header->SetCentrality(0);
-  }
-  
-  //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()};
-  Float_t diamcov[3];
-  event->GetDiamondCovXY(diamcov);
-  header->SetDiamond(diamxy,diamcov);
-  header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
-  
-}
-
-//_____________________________________________
-void AliAnalysisTaskCaloFilter::FillAODTracks()
-{
-  // AOD track copy
-  
-  if(!fFillTracks) return;
-  
-  AliVEvent* event = InputEvent();
-  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);  
-  //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
-
-  // Copy from AODs
-  if(aodevent)
+  // AODs
+  if(!fESDEvent) 
   {
-    //TClonesArray* inTracks = aodevent  ->GetTracks();
-    TClonesArray* ouTracks = AODEvent()->GetTracks();
-               //new (ouTracks) TClonesArray(*inTracks);
-        
-    for (Int_t nTrack = 0; nTrack < aodevent->GetNumberOfTracks(); ++nTrack) 
+    // 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 ) 
     {
-      AliAODTrack *track = aodevent->GetTrack(nTrack);
-
-      new((*ouTracks)[nTrack]) AliAODTrack(*track);
-    } 
+      if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject \n");
+      
+      return kFALSE ;
+    }
     
-    return;
+    return kTRUE;
   }
   
-}  
-
-//_______________________________________________
-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.;
+  // ESDs
+  if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() > 0)
+  {
+    return kTRUE;
+  }
   
-  // Copy from AODs
-  if(aodevent)
+  if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() < 1) 
   {
-    TClonesArray* inVertices = aodevent  ->GetVertices();
-    TClonesArray* ouVertices = AODEvent()->GetVertices();
-
-               //new (ouVertices) TClonesArray(*inVertices);
-    
-    for (Int_t nVertices = 0; nVertices < inVertices->GetEntriesFast(); ++nVertices) 
+    // SPD vertex
+    if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() > 0) 
     {
-      AliAODVertex *vertex = (AliAODVertex*)inVertices->At(nVertices);
+      //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+      return kTRUE;
       
-      new((*ouVertices)[nVertices]) AliAODVertex(*vertex);
-    } 
-    
-    return;
+    }
+    if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() < 1) 
+    {
+      //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+      if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject \n");
+      return kFALSE;
+    }
   }
   
-  if(!esdevent) return;
+  if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject \n");
   
-  // Copy from ESDs
-  
-  // Access to the AOD container of vertices
-  Int_t jVertices=0;
-  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();
-  
-  AliAODVertex * primary = new(vertices[jVertices++])
-  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
-  primary->SetName(event->GetPrimaryVertex()->GetName());
-  primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
+  return kFALSE;
   
 }
 
@@ -653,65 +419,594 @@ void AliAnalysisTaskCaloFilter::CorrectionsInEMCAL()
   } // corrections in EMCAL
 }
 
-//_______________________________________________________
-Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex()
+//________________________________________________
+void AliAnalysisTaskCaloFilter::FillAODCaloCells()
+{  
+  // Fill EMCAL/PHOS cell info
+  
+  // EMCAL
+  if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && fEvent->GetEMCALCells()) 
+  { // protection against missing ESD information
+    AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
+    Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
+    
+    AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
+    aodEMcells.CreateContainer(nEMcell);
+    aodEMcells.SetType(AliVCaloCells::kEMCALCell);
+    Double_t calibFactor = 1.;   
+    for (Int_t iCell = 0; iCell < nEMcell; iCell++) 
+    { 
+      Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
+      fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); 
+      fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);      
+      
+      if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn())
+      { 
+        calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
+      }
+      
+      if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
+      { //Channel is not declared as bad
+        aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
+                           eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
+        //printf("GOOD channel\n");
+      }
+      else 
+      {
+        aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
+        //printf("BAD channel\n");
+      }
+    }
+    aodEMcells.Sort();
+  }
+  
+  // PHOS 
+  if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && fEvent->GetPHOSCells()) 
+  { // protection against missing ESD information
+    AliVCaloCells &eventPHcells = *(fEvent->GetPHOSCells());
+    Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
+    
+    AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
+    aodPHcells.CreateContainer(nPHcell);
+    aodPHcells.SetType(AliVCaloCells::kPHOSCell);
+    
+    for (Int_t iCell = 0; iCell < nPHcell; iCell++) 
+    {      
+      aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell),
+                         eventPHcells.GetTime(iCell),eventPHcells.GetMCLabel(iCell),eventPHcells.GetEFraction(iCell));
+    }
+    
+    aodPHcells.Sort();
+  }
+}
+
+
+//___________________________________________________
+void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
 {
-  //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
-  //It only works for ESDs
+  // Fill the AOD with caloclusters
   
-  AliESDEvent * event = dynamic_cast<AliESDEvent*> (InputEvent());
+  // Access to the AOD container of clusters
   
-  // AODs
-  if(!event) 
+  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
+  Int_t jClusters=0;
+  Float_t  posF[3]  ;
+  
+  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
+  for (Int_t iClust=0; iClust<nCaloClus; ++iClust) 
   {
-    // 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 ;
+    AliVCluster * cluster = fEvent->GetCaloCluster(iClust);
     
-    return kTRUE;
+    //Check which calorimeter information we want to keep.
+    
+    if(fCaloFilter!=kBoth)
+    {
+      if     (fCaloFilter==kPHOS  && cluster->IsEMCAL()) continue ;
+      else if(fCaloFilter==kEMCAL && cluster->IsPHOS())  continue ;
+    }  
+    
+    // Get original residuals, in case of previous recalculation, reset them 
+    Float_t dR = cluster->GetTrackDx();
+    Float_t dZ = cluster->GetTrackDz();
+    
+    if(DebugLevel() > 2)
+      printf("Original residuals : dZ %f, dR %f\n ",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());
+      
+      if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;       
+      
+      if(fEMCALRecoUtils->IsExoticCluster(cluster, InputEvent()->GetEMCALCells(),InputEvent()->GetBunchCrossNumber())) continue;       
+      
+      fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
+      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);
+    }
+    
+    //--------------------------------------------------------------
+    
+    //Now fill AODs
+    
+    Int_t id       = cluster->GetID();
+    Float_t energy = cluster->E();
+    cluster->GetPosition(posF);
+    
+    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
+    AliAODCaloCluster(id,
+                      cluster->GetNLabels(),
+                      cluster->GetLabels(),
+                      energy,
+                      posF,
+                      NULL,
+                      cluster->GetType());
+    
+    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
+                                cluster->GetDispersion(),
+                                cluster->GetM20(), cluster->GetM02(),
+                                -1,  
+                                cluster->GetNExMax(),cluster->GetTOF()) ;
+    
+    caloCluster->SetPIDFromESD(cluster->GetPID());
+    caloCluster->SetNCells(cluster->GetNCells());
+    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
+    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]);
+    }    
+    
+    //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) 
+    { //Default value in PHOS 999, in EMCAL 1024, why?
+      caloCluster->AddTrackMatched(new AliAODTrack); 
+    }
+    // TO DO, in case Tracks available, think how to put the matched track in AOD
   }
   
-  // ESDs
-  if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
+  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots  
+  
+}
+
+//__________________________________________________
+void AliAnalysisTaskCaloFilter::FillAODCaloTrigger()
+{
+  // AOD CaloTrigger copy
+  
+  AliAODCaloTrigger* triggerEM = AODEvent()->GetCaloTrigger("EMCAL");
+  AliAODCaloTrigger* triggerPH = AODEvent()->GetCaloTrigger("PHOS");
+  
+  // Copy from AODs
+  if(fAODEvent)
   {
-    return kTRUE;
+    if(fCaloFilter==kBoth || fCaloFilter==kPHOS)  *triggerPH = *(fAODEvent->GetCaloTrigger("PHOS"));
+    
+    if(fCaloFilter==kBoth || fCaloFilter==kEMCAL) *triggerEM = *(fAODEvent->GetCaloTrigger("EMCAL"));
+    
+    return;
   }
   
-  if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) 
+}  
+
+//______________________________________________
+void AliAnalysisTaskCaloFilter::FillAODHeader()
+{
+  // AOD header copy
+  
+  AliAODHeader* header = AODEvent()->GetHeader();
+  
+  // Copy from AODs
+  if(fAODEvent)
   {
-    // SPD vertex
-    if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) 
+    *header = *(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(fEvent->GetBunchCrossNumber());
+  header->SetOrbitNumber(fEvent->GetOrbitNumber());
+  header->SetPeriodNumber(fEvent->GetPeriodNumber());
+  header->SetEventType(fEvent->GetEventType());
+  
+  //Centrality
+  if(fEvent->GetCentrality())
+  {
+    header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
+  }
+  else
+  {
+    header->SetCentrality(0);
+  }
+  
+  //Trigger  
+  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
+  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];
+  fEvent->GetDiamondCovXY(diamcov);
+  header->SetDiamond(diamxy,diamcov);
+  header->SetDiamondZ(fESDEvent->GetDiamondZ(),fESDEvent->GetSigma2DiamondZ());
+  
+}
+
+//_____________________________________________
+void AliAnalysisTaskCaloFilter::FillAODTracks()
+{
+  // AOD track copy
+  
+  if(!fFillTracks) return;
+  
+  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(fAODEvent)
+  {
+    //TClonesArray* inTracks = fAODEvent  ->GetTracks();
+    TClonesArray* ouTracks = AODEvent()->GetTracks();
+               //new (ouTracks) TClonesArray(*inTracks);
+    
+    //printf("N tracks %d\n",fAODEvent->GetNumberOfTracks());
+    Int_t nCopyTrack = 0;
+    for (Int_t nTrack = 0; nTrack < fAODEvent->GetNumberOfTracks(); ++nTrack) 
     {
-      //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
-      return kTRUE;
+      AliAODTrack *track = fAODEvent->GetTrack(nTrack);
       
-    }
-    if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) 
+      // 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(),
+                                                            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) 
     {
-      //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
-      return kFALSE;
+      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
+  
+  // 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(fAODEvent)
+  {
+    TClonesArray* inVertices = fAODEvent  ->GetVertices();
+    TClonesArray* ouVertices = AODEvent()->GetVertices();
+    
+               //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;
   }
-  return kFALSE;
+  
+  if(!fESDEvent) return;
+  
+  // Copy from ESDs
+  
+  // Access to the AOD container of vertices
+  Int_t jVertices=0;
+  TClonesArray &vertices = *(AODEvent()->GetVertices());
+  
+  // Add primary vertex. The primary tracks will be defined
+  // 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(fEvent->GetPrimaryVertex()->GetName());
+  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
   
 }
 
+//____________________________________
+void AliAnalysisTaskCaloFilter::Init()
+{
+  //Init analysis with configuration macro if available
+  
+  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;
+    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;
+    
+    for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
+  }
+} 
+
 //_________________________________________
 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 8f999f8a14e204d33dd5da42901afa264e38fc66..f75055efa81b5199d8615fa32841ee6fa7c3c8b4 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 c5758866ba0dbaea5cb608abf61471418ed631fd..a1c122f44c527e74faa7c4475ae96c5677bf3de1 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();