[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskCaloFilter.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id: AliAnalysisTaskCaloFilter.cxx $ */

//////////////////////////////////////////////////////////
// Filter the ESDCaloClusters and ESDCaloCells of EMCAL,
// PHOS or both, creating the corresponing AODCaloClusters
// and AODCaloCells.
// Keep also the AODHeader information and the vertex.
// Needed for calorimeter calibration.
// Copy of AliAnalysisTaskESDfilter.
// Author: Gustavo Conesa Balbastre (INFN - Frascati)
//////////////////////////////////////////////////////////

#include "AliAnalysisTaskCaloFilter.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliLog.h"
#include "AliVCluster.h"
#include "AliVCaloCells.h"
#include "AliEMCALRecoUtils.h"
#include "AliEMCALGeometry.h"
#include "AliVEventHandler.h"
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliESDtrackCuts.h"

ClassImp(AliAnalysisTaskCaloFilter)
  
////////////////////////////////////////////////////////////////////////

AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
  AliAnalysisTaskSE(), //fCuts(0x0),
  fCaloFilter(0), fCorrect(kFALSE), 
  fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEAR"), 
  fEMCALRecoUtils(new AliEMCALRecoUtils),
  fESDtrackCuts(0), fTrackMultEtaCut(0.9)
{
  // Default constructor
  fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();

}

//__________________________________________________
AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
  AliAnalysisTaskSE(name), //fCuts(0x0),
  fCaloFilter(0), fCorrect(kFALSE),
  fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEAR"), 
  fEMCALRecoUtils(new AliEMCALRecoUtils),
  fESDtrackCuts(0), fTrackMultEtaCut(0.9)
{
  // Constructor
  fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();

}

//__________________________________________________
AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
{
  //Destructor.
	
  if(fEMCALGeo)       delete fEMCALGeo;	
  if(fEMCALRecoUtils) delete fEMCALRecoUtils;
  if(fESDtrackCuts)   delete fESDtrackCuts;
}

//__________________________________________________
void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
{
  // Init geometry 
	
  fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;	
  
}  

//__________________________________________________
void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event
  //
  
  if (fDebug > 0)  
    printf("CaloFilter: Analysing event # %d\n", (Int_t)Entry());
  
  // Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
  
  AliVEvent* event = InputEvent();
  if(!event) {
    printf("AliAnalysisTaskCaloFilter::CreateAODFromESD() - This event does not contain Input?");
    return;
  }

  //Magic line to write events to file
  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
    
  Bool_t bAOD = kFALSE;
  if(!strcmp(event->GetName(),"AliAODEvent")) bAOD=kTRUE;
  Bool_t bESD = kFALSE;
  if(!strcmp(event->GetName(),"AliESDEvent")) bESD=kTRUE;
  
  //Get track multiplicity
  Int_t trackMult = 0;
  if(bESD){
    Int_t nTracks   = InputEvent()->GetNumberOfTracks() ;
    for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
      AliVTrack * track = (AliVTrack*)InputEvent()->GetTrack(itrack) ; // retrieve track from esd
      if(!fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
      //Count the tracks in eta < 0.9
      if(TMath::Abs(track->Eta())< fTrackMultEtaCut) trackMult++;
    }    
  }
  
  // set arrays and pointers
  Float_t posF[3];
  Double_t pos[3];
  
  Double_t covVtx[6];
  
  for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
      
  AliAODHeader* header = AODEvent()->GetHeader();
  
  header->SetRunNumber(event->GetRunNumber());
  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
  if(bESD)
    header->SetFiredTriggerClasses(((AliESDEvent*)event)->GetFiredTriggerClasses());
  header->SetTriggerMask(event->GetTriggerMask()); 
  header->SetTriggerCluster(event->GetTriggerCluster());
  
  header->SetBunchCrossNumber(event->GetBunchCrossNumber());
  header->SetOrbitNumber(event->GetOrbitNumber());
  header->SetPeriodNumber(event->GetPeriodNumber());
  header->SetEventType(event->GetEventType());
  header->SetMuonMagFieldScale(-999.); // FIXME
  //printf("Track Multiplicity for eta < %f: %d \n",fTrackMultEtaCut,trackMult);
  header->SetCentrality((Double_t)trackMult);        // FIXME
  //printf("Centrality %f\n",header->GetCentrality());
  
  header->SetTriggerMask(event->GetTriggerMask()); 
  header->SetTriggerCluster(event->GetTriggerCluster());
  header->SetMagneticField(event->GetMagneticField());
  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);
  if(bESD){
    header->SetDiamondZ(((AliESDEvent*)event)->GetDiamondZ(),((AliESDEvent*)event)->GetSigma2DiamondZ());
  }
  //
  //
  Int_t nVertices = 1 ;/* = prim. vtx*/;
  Int_t nCaloClus = event->GetNumberOfCaloClusters();
  
  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
  
  // Access to the AOD container of vertices
  TClonesArray &vertices = *(AODEvent()->GetVertices());
  Int_t jVertices=0;
  
  // Add primary vertex. The primary tracks will be defined
  // after the loops on the composite objects (V0, cascades, kinks)
  event->GetPrimaryVertex()->GetXYZ(pos);
  Float_t chi = 0;
  if      (bESD){
    ((AliESDEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
    chi = ((AliESDEvent*)event)->GetPrimaryVertex()->GetChi2toNDF();
  }
  else if (bAOD){
    ((AliAODEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
    chi = ((AliAODEvent*)event)->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
  }
  
  AliAODVertex * primary = new(vertices[jVertices++])
    AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
  primary->SetName(event->GetPrimaryVertex()->GetName());
  primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
  
  // Access to the AOD container of clusters
  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
  Int_t jClusters=0;
  
  for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
    
    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 ;
    }  
    
    //--------------------------------------------------------------
    //If EMCAL, and requested, correct energy, position ...
    if(cluster->IsEMCAL() && fCorrect){
      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(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
//        printf("Finally reject\n");
//        continue;
//      }
      if(DebugLevel() > 2)
      { 
        printf("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClust,cluster->E(),
               cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
        cluster->GetPosition(position);
        printf("Filter, before  : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
      }
            
      if(fEMCALRecoUtils->IsRecalibrationOn())	{
        fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
        fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, event->GetEMCALCells(),cluster);
        fEMCALRecoUtils->RecalculateClusterPID(cluster);

      }
      cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
      
      fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);

      if(DebugLevel() > 2)
      { 
        printf("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",cluster->GetID(),cluster->E(),
               cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
        cluster->GetPosition(position);
        printf("Filter, after   : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
      }    
      
    }
    //--------------------------------------------------------------

    //Now fill AODs
    Int_t id       = cluster->GetID();
    Float_t energy = cluster->E();
    cluster->GetPosition(posF);
    
    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
      AliAODCaloCluster(id,
			0,
			0x0,
			energy,
			posF,
			NULL,
			cluster->GetType());
    
    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
				cluster->GetDispersion(),
				cluster->GetM20(), cluster->GetM02(),
				cluster->GetEmcCpvDistance(),  
				cluster->GetNExMax(),cluster->GetTOF()) ;
    
    caloCluster->SetPIDFromESD(cluster->GetPID());
    caloCluster->SetNCells(cluster->GetNCells());
    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
    
    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(bESD){
      TArrayI* matchedT = 	((AliESDCaloCluster*)cluster)->GetTracksMatched();
      if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {	
        for (Int_t im = 0; im < matchedT->GetSize(); im++) {
          Int_t iESDtrack = matchedT->At(im);;
          if ((AliVTrack*)InputEvent()->GetTrack(iESDtrack) != 0) {
            caloCluster->AddTrackMatched((AliVTrack*)InputEvent()->GetTrack(iESDtrack));
          }
        }
      }// There is at least a match with a track
    }
  } 
  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters	 
  // end of loop on calo clusters
  
  // fill EMCAL cell info
  if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && event->GetEMCALCells()) { // protection against missing ESD information
    AliVCaloCells &eventEMcells = *(event->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);
        //printf("GOOD channel\n");
      }
      else {
        aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
        //printf("BAD channel\n");

      }
    }
    aodEMcells.Sort();
  }
  
  // fill PHOS cell info
  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));
    }
    aodPHcells.Sort();
  }
  
  
  return;
}

//_____________________________________________________
void AliAnalysisTaskCaloFilter::PrintInfo(){

  //Print settings

  printf("TASK: AnalysisCaloFilter \n");
  printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
  printf("\t Calorimeter Filtering Option     ? %d\n",fCaloFilter);

}

//_____________________________________________________
//void AliAnalysisTaskCaloFilter::LocalInit()
//{
//	// Local Initialization
//	
//	// Create cuts/param objects and publish to slot
//	const Int_t buffersize = 255;
//	char onePar[buffersize] ;
//	fCuts = new TList();
//  
//	snprintf(onePar,buffersize, "Calorimeter Filtering Option %d", fCaloFilter) ;
//	fCuts->Add(new TObjString(onePar));
//	snprintf(onePar,buffersize, "Not only filter but correct? %d cells;", fCorrect) ;
//	fCuts->Add(new TObjString(onePar));
//  
//	// Post Data
//	PostData(2, fCuts);
//	
//}


//__________________________________________________
void AliAnalysisTaskCaloFilter::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
    if (fDebug > 1) printf("AnalysisCaloFilter: Terminate() \n");
}
Commit	Line	Data
7a4cf423	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id: AliAnalysisTaskCaloFilter.cxx $ */
	17
	18	//////////////////////////////////////////////////////////
	19	// Filter the ESDCaloClusters and ESDCaloCells of EMCAL,
	20	// PHOS or both, creating the corresponing AODCaloClusters
	21	// and AODCaloCells.
	22	// Keep also the AODHeader information and the vertex.
	23	// Needed for calorimeter calibration.
	24	// Copy of AliAnalysisTaskESDfilter.
	25	// Author: Gustavo Conesa Balbastre (INFN - Frascati)
	26	//////////////////////////////////////////////////////////
	27
	28	#include "AliAnalysisTaskCaloFilter.h"
	29	#include "AliESDEvent.h"
	30	#include "AliAODEvent.h"
7a4cf423	31	#include "AliLog.h"
247abff4	32	#include "AliVCluster.h"
	33	#include "AliVCaloCells.h"
	34	#include "AliEMCALRecoUtils.h"
	35	#include "AliEMCALGeometry.h"
	36	#include "AliVEventHandler.h"
	37	#include "AliAnalysisManager.h"
	38	#include "AliInputEventHandler.h"
3a58eee6	39	#include "AliESDtrackCuts.h"
7a4cf423	40
	41	ClassImp(AliAnalysisTaskCaloFilter)
	42
	43	////////////////////////////////////////////////////////////////////////
	44
	45	AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
247abff4	46	AliAnalysisTaskSE(), //fCuts(0x0),
	47	fCaloFilter(0), fCorrect(kFALSE),
	48	fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEAR"),
3a58eee6	49	fEMCALRecoUtils(new AliEMCALRecoUtils),
3a58eee6	50	fESDtrackCuts(0), fTrackMultEtaCut(0.9)
7a4cf423	51	{
7a4cf423	52	// Default constructor
3a58eee6	53	fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
3a58eee6	54
7a4cf423	55	}
	56
	57	//__________________________________________________
	58	AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
247abff4	59	AliAnalysisTaskSE(name), //fCuts(0x0),
	60	fCaloFilter(0), fCorrect(kFALSE),
	61	fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEAR"),
3a58eee6	62	fEMCALRecoUtils(new AliEMCALRecoUtils),
3a58eee6	63	fESDtrackCuts(0), fTrackMultEtaCut(0.9)
7a4cf423	64	{
7a4cf423	65	// Constructor
3a58eee6	66	fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
3a58eee6	67
7a4cf423	68	}
	69
	70	//__________________________________________________
247abff4	71	AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
7a4cf423	72	{
247abff4	73	//Destructor.
	74
	75	if(fEMCALGeo) delete fEMCALGeo;
3a58eee6	76	if(fEMCALRecoUtils) delete fEMCALRecoUtils;
3a58eee6	77	if(fESDtrackCuts) delete fESDtrackCuts;
247abff4	78	}
	79
	80	//__________________________________________________
	81	void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
	82	{
	83	// Init geometry
	84
	85	fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
7a4cf423	86
247abff4	87	}
	88
	89	//__________________________________________________
	90	void AliAnalysisTaskCaloFilter::UserExec(Option_t /option*/)
	91	{
	92	// Execute analysis for current event
7a4cf423	93	//
7a4cf423	94
247abff4	95	if (fDebug > 0)
247abff4	96	printf("CaloFilter: Analysing event # %d\n", (Int_t)Entry());
7a4cf423	97
247abff4	98	// Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
7a4cf423	99
247abff4	100	AliVEvent* event = InputEvent();
247abff4	101	if(!event) {
247abff4	102	printf("AliAnalysisTaskCaloFilter::CreateAODFromESD() - This event does not contain Input?");
2dfb1428	103	return;
2dfb1428	104	}
fbad5435	105
	106	//Magic line to write events to file
	107	AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
	108
fbad5435	109	Bool_t bAOD = kFALSE;
	110	if(!strcmp(event->GetName(),"AliAODEvent")) bAOD=kTRUE;
	111	Bool_t bESD = kFALSE;
	112	if(!strcmp(event->GetName(),"AliESDEvent")) bESD=kTRUE;
3a58eee6	113
	114	//Get track multiplicity
	115	Int_t trackMult = 0;
	116	if(bESD){
	117	Int_t nTracks = InputEvent()->GetNumberOfTracks() ;
	118	for (Int_t itrack = 0; itrack < nTracks; itrack++) {////////////// track loop
	119	AliVTrack * track = (AliVTrack*)InputEvent()->GetTrack(itrack) ; // retrieve track from esd
	120	if(!fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
	121	//Count the tracks in eta < 0.9
	122	if(TMath::Abs(track->Eta())< fTrackMultEtaCut) trackMult++;
	123	}
	124	}
	125
7a4cf423	126	// set arrays and pointers
	127	Float_t posF[3];
	128	Double_t pos[3];
	129
	130	Double_t covVtx[6];
	131
	132	for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
247abff4	133
7a4cf423	134	AliAODHeader* header = AODEvent()->GetHeader();
247abff4	135
	136	header->SetRunNumber(event->GetRunNumber());
	137	header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
	138	if(bESD)
	139	header->SetFiredTriggerClasses(((AliESDEvent*)event)->GetFiredTriggerClasses());
	140	header->SetTriggerMask(event->GetTriggerMask());
	141	header->SetTriggerCluster(event->GetTriggerCluster());
	142
	143	header->SetBunchCrossNumber(event->GetBunchCrossNumber());
	144	header->SetOrbitNumber(event->GetOrbitNumber());
	145	header->SetPeriodNumber(event->GetPeriodNumber());
	146	header->SetEventType(event->GetEventType());
7a4cf423	147	header->SetMuonMagFieldScale(-999.); // FIXME
3a58eee6	148	//printf("Track Multiplicity for eta < %f: %d \n",fTrackMultEtaCut,trackMult);
	149	header->SetCentrality((Double_t)trackMult); // FIXME
	150	//printf("Centrality %f\n",header->GetCentrality());
7a4cf423	151
247abff4	152	header->SetTriggerMask(event->GetTriggerMask());
	153	header->SetTriggerCluster(event->GetTriggerCluster());
	154	header->SetMagneticField(event->GetMagneticField());
	155	header->SetZDCN1Energy(event->GetZDCN1Energy());
	156	header->SetZDCP1Energy(event->GetZDCP1Energy());
	157	header->SetZDCN2Energy(event->GetZDCN2Energy());
	158	header->SetZDCP2Energy(event->GetZDCP2Energy());
	159	header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
	160	Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
	161	Float_t diamcov[3]; event->GetDiamondCovXY(diamcov);
7a4cf423	162	header->SetDiamond(diamxy,diamcov);
247abff4	163	if(bESD){
	164	header->SetDiamondZ(((AliESDEvent)event)->GetDiamondZ(),((AliESDEvent)event)->GetSigma2DiamondZ());
	165	}
7a4cf423	166	//
	167	//
	168	Int_t nVertices = 1 ;/* = prim. vtx*/;
247abff4	169	Int_t nCaloClus = event->GetNumberOfCaloClusters();
7a4cf423	170
	171	AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
	172
	173	// Access to the AOD container of vertices
	174	TClonesArray &vertices = *(AODEvent()->GetVertices());
	175	Int_t jVertices=0;
	176
	177	// Add primary vertex. The primary tracks will be defined
	178	// after the loops on the composite objects (V0, cascades, kinks)
247abff4	179	event->GetPrimaryVertex()->GetXYZ(pos);
	180	Float_t chi = 0;
	181	if (bESD){
	182	((AliESDEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
	183	chi = ((AliESDEvent*)event)->GetPrimaryVertex()->GetChi2toNDF();
	184	}
	185	else if (bAOD){
	186	((AliAODEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
	187	chi = ((AliAODEvent*)event)->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
	188	}
7a4cf423	189
7a4cf423	190	AliAODVertex * primary = new(vertices[jVertices++])
247abff4	191	AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
	192	primary->SetName(event->GetPrimaryVertex()->GetName());
	193	primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
7a4cf423	194
	195	// Access to the AOD container of clusters
	196	TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
	197	Int_t jClusters=0;
	198
	199	for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
	200
247abff4	201	AliVCluster * cluster = event->GetCaloCluster(iClust);
7a4cf423	202
7a4cf423	203	//Check which calorimeter information we want to keep.
7a4cf423	204
247abff4	205	if(fCaloFilter!=kBoth){
	206	if (fCaloFilter==kPHOS && cluster->IsEMCAL()) continue ;
	207	else if(fCaloFilter==kEMCAL && cluster->IsPHOS()) continue ;
	208	}
	209
	210	//--------------------------------------------------------------
	211	//If EMCAL, and requested, correct energy, position ...
247abff4	212	if(cluster->IsEMCAL() && fCorrect){
	213	Float_t position[]={0,0,0};
	214	if(DebugLevel() > 2)
	215	printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
	216	if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
5ef94e1b	217	// if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
	218	// printf("Finally reject\n");
	219	// continue;
	220	// }
247abff4	221	if(DebugLevel() > 2)
	222	{
	223	printf("Filter, before : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClust,cluster->E(),
	224	cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
	225	cluster->GetPosition(position);
	226	printf("Filter, before : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
	227	}
	228
5ef94e1b	229	if(fEMCALRecoUtils->IsRecalibrationOn()) {
	230	fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
	231	fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, event->GetEMCALCells(),cluster);
	232	fEMCALRecoUtils->RecalculateClusterPID(cluster);
	233
	234	}
247abff4	235	cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
247abff4	236
5ef94e1b	237	fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);
5ef94e1b	238
247abff4	239	if(DebugLevel() > 2)
	240	{
	241	printf("Filter, after : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",cluster->GetID(),cluster->E(),
	242	cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
	243	cluster->GetPosition(position);
	244	printf("Filter, after : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
	245	}
	246
	247	}
	248	//--------------------------------------------------------------
	249
	250	//Now fill AODs
7a4cf423	251	Int_t id = cluster->GetID();
	252	Float_t energy = cluster->E();
	253	cluster->GetPosition(posF);
7a4cf423	254
	255	AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
	256	AliAODCaloCluster(id,
	257	0,
	258	0x0,
	259	energy,
	260	posF,
	261	NULL,
c8fe2783	262	cluster->GetType());
7a4cf423	263
7a4cf423	264	caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
c8fe2783	265	cluster->GetDispersion(),
7a4cf423	266	cluster->GetM20(), cluster->GetM02(),
	267	cluster->GetEmcCpvDistance(),
	268	cluster->GetNExMax(),cluster->GetTOF()) ;
	269
c8fe2783	270	caloCluster->SetPIDFromESD(cluster->GetPID());
7a4cf423	271	caloCluster->SetNCells(cluster->GetNCells());
	272	caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
	273	caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
	274
247abff4	275	if(DebugLevel() > 2)
	276	{
	277	printf("Filter, aod : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
	278	caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
	279	caloCluster->GetPosition(posF);
	280	printf("Filter, aod : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
	281	}
	282
3a58eee6	283	//Matched tracks, just to know if there was any match, the track pointer is useless.
	284	if(bESD){
	285	TArrayI* matchedT = ((AliESDCaloCluster*)cluster)->GetTracksMatched();
	286	if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
	287	for (Int_t im = 0; im < matchedT->GetSize(); im++) {
	288	Int_t iESDtrack = matchedT->At(im);;
	289	if ((AliVTrack*)InputEvent()->GetTrack(iESDtrack) != 0) {
	290	caloCluster->AddTrackMatched((AliVTrack*)InputEvent()->GetTrack(iESDtrack));
	291	}
	292	}
	293	}// There is at least a match with a track
	294	}
7a4cf423	295	}
	296	caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
	297	// end of loop on calo clusters
	298
	299	// fill EMCAL cell info
247abff4	300	if ((fCaloFilter==kBoth \|\| fCaloFilter==kEMCAL) && event->GetEMCALCells()) { // protection against missing ESD information
	301	AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
	302	Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
7a4cf423	303
	304	AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
	305	aodEMcells.CreateContainer(nEMcell);
c8fe2783	306	aodEMcells.SetType(AliVCaloCells::kEMCALCell);
247abff4	307	Double_t calibFactor = 1.;
247abff4	308	for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
5ef94e1b	309	Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
	310	fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
	311	fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
	312
247abff4	313	if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn()){
247abff4	314	calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
247abff4	315	}
5ef94e1b	316
	317	if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
	318	aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
	319	//printf("GOOD channel\n");
	320	}
	321	else {
	322	aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
	323	//printf("BAD channel\n");
	324
	325	}
7a4cf423	326	}
	327	aodEMcells.Sort();
	328	}
	329
	330	// fill PHOS cell info
247abff4	331	if ((fCaloFilter==kBoth \|\| fCaloFilter==kPHOS) && event->GetPHOSCells()) { // protection against missing ESD information
	332	AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
	333	Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
7a4cf423	334
	335	AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
	336	aodPHcells.CreateContainer(nPHcell);
c8fe2783	337	aodPHcells.SetType(AliVCaloCells::kPHOSCell);
7a4cf423	338	for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
247abff4	339	aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell));
7a4cf423	340	}
	341	aodPHcells.Sort();
	342	}
	343
	344
	345	return;
	346	}
	347
5ef94e1b	348	//_____________________________________________________
	349	void AliAnalysisTaskCaloFilter::PrintInfo(){
	350
	351	//Print settings
	352
	353	printf("TASK: AnalysisCaloFilter \n");
	354	printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
	355	printf("\t Calorimeter Filtering Option ? %d\n",fCaloFilter);
	356
	357	}
	358
247abff4	359	//_____________________________________________________
	360	//void AliAnalysisTaskCaloFilter::LocalInit()
	361	//{
	362	// // Local Initialization
	363	//
	364	// // Create cuts/param objects and publish to slot
	365	// const Int_t buffersize = 255;
	366	// char onePar[buffersize] ;
	367	// fCuts = new TList();
	368	//
	369	// snprintf(onePar,buffersize, "Calorimeter Filtering Option %d", fCaloFilter) ;
	370	// fCuts->Add(new TObjString(onePar));
	371	// snprintf(onePar,buffersize, "Not only filter but correct? %d cells;", fCorrect) ;
	372	// fCuts->Add(new TObjString(onePar));
	373	//
	374	// // Post Data
	375	// PostData(2, fCuts);
	376	//
	377	//}
7a4cf423	378
7a4cf423	379
	380	//__________________________________________________
	381	void AliAnalysisTaskCaloFilter::Terminate(Option_t /option*/)
	382	{
	383	// Terminate analysis
	384	//
	385	if (fDebug > 1) printf("AnalysisCaloFilter: Terminate() \n");
	386	}
	387