[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;
  
  //Do Corrections in EMCAL 
  //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) ) {
    //Cluster Loop
    for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
      
      AliVCluster * cluster = event->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(!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]);
      }    
      
    }
    //Recalculate track-matching
    fEMCALRecoUtils->FindMatches(event);
    
  } // corrections in EMCAL
  
  //Now loop on clusters to fill AODs
  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 ;
    }  
    
    //Temporary trick, FIXME
    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->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
      //        printf("Finally reject\n");
      //        continue;
      //      }
      
      fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
      if(DebugLevel() > 2)
        printf("Corrected 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,
			0,
			0x0,
			energy,
			posF,
			NULL,
			cluster->GetType());
    
    caloCluster->SetChi2(dZ);//Temporary trick, FIXME
    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
				cluster->GetDispersion(),
				cluster->GetM20(), cluster->GetM02(),
				dR,  //Temporary trick, FIXME
				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.
    //Temporary trick, FIXME
    if(bESD){
      if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990) { //Default value in PHOS 999, in EMCAL 1024, why?
        if(DebugLevel() > 2) 
          printf("*** Cluster Track-Matched *** dR %f, dZ %f\n",caloCluster->GetEmcCpvDistance(),caloCluster->Chi2());
        caloCluster->AddTrackMatched(0x0); 
      }// fill the array with one entry to signal a possible match
      //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
f2ccb5b8	199	//Do Corrections in EMCAL
	200	//If EMCAL, and requested, correct energy, position ...
	201	//Need to do this in a separate loop before filling the ESDs because of the track matching recalculations
	202	if(fCorrect && (fCaloFilter==kEMCAL \|\| fCaloFilter==kBoth) ) {
	203	//Cluster Loop
	204	for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
	205
	206	AliVCluster * cluster = event->GetCaloCluster(iClust);
	207	if(cluster->IsPHOS()) continue ;
	208
247abff4	209	Float_t position[]={0,0,0};
	210	if(DebugLevel() > 2)
	211	printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
	212	if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
f2ccb5b8	213	// if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
	214	// printf("Finally reject\n");
	215	// continue;
	216	// }
247abff4	217	if(DebugLevel() > 2)
	218	{
	219	printf("Filter, before : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClust,cluster->E(),
	220	cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
	221	cluster->GetPosition(position);
	222	printf("Filter, before : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
	223	}
f2ccb5b8	224
5ef94e1b	225	if(fEMCALRecoUtils->IsRecalibrationOn()) {
	226	fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
	227	fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, event->GetEMCALCells(),cluster);
	228	fEMCALRecoUtils->RecalculateClusterPID(cluster);
5ef94e1b	229	}
247abff4	230	cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
247abff4	231
5ef94e1b	232	fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);
f2ccb5b8	233
247abff4	234	if(DebugLevel() > 2)
	235	{
	236	printf("Filter, after : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",cluster->GetID(),cluster->E(),
	237	cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
	238	cluster->GetPosition(position);
	239	printf("Filter, after : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
	240	}
	241
f2ccb5b8	242	}
	243	//Recalculate track-matching
	244	fEMCALRecoUtils->FindMatches(event);
	245
	246	} // corrections in EMCAL
	247
	248	//Now loop on clusters to fill AODs
	249	for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
	250
	251	AliVCluster * cluster = event->GetCaloCluster(iClust);
	252
	253	//Check which calorimeter information we want to keep.
	254
	255	if(fCaloFilter!=kBoth){
	256	if (fCaloFilter==kPHOS && cluster->IsEMCAL()) continue ;
	257	else if(fCaloFilter==kEMCAL && cluster->IsPHOS()) continue ;
	258	}
	259
	260	//Temporary trick, FIXME
	261	Float_t dR = cluster->GetTrackDx();
	262	Float_t dZ = cluster->GetTrackDz();
	263	if(DebugLevel() > 2)
	264	printf("Original residuals : dZ %f, dR %f\n ",dZ, dR);
	265	//--------------------------------------------------------------
	266	//If EMCAL and corrections done, get the new matching parameters, do not copy noisy clusters
	267	if(cluster->IsEMCAL() && fCorrect){
	268	if(DebugLevel() > 2)
	269	printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
	270	if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
	271	// if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
	272	// printf("Finally reject\n");
	273	// continue;
	274	// }
	275
	276	fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
	277	if(DebugLevel() > 2)
	278	printf("Corrected Residuals : dZ %f, dR %f\n ",dZ, dR);
	279
247abff4	280	}
	281	//--------------------------------------------------------------
	282
	283	//Now fill AODs
f2ccb5b8	284
7a4cf423	285	Int_t id = cluster->GetID();
	286	Float_t energy = cluster->E();
	287	cluster->GetPosition(posF);
7a4cf423	288
	289	AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
	290	AliAODCaloCluster(id,
	291	0,
	292	0x0,
	293	energy,
	294	posF,
	295	NULL,
c8fe2783	296	cluster->GetType());
7a4cf423	297
f2ccb5b8	298	caloCluster->SetChi2(dZ);//Temporary trick, FIXME
7a4cf423	299	caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
c8fe2783	300	cluster->GetDispersion(),
7a4cf423	301	cluster->GetM20(), cluster->GetM02(),
f2ccb5b8	302	dR, //Temporary trick, FIXME
7a4cf423	303	cluster->GetNExMax(),cluster->GetTOF()) ;
7a4cf423	304
c8fe2783	305	caloCluster->SetPIDFromESD(cluster->GetPID());
7a4cf423	306	caloCluster->SetNCells(cluster->GetNCells());
	307	caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
	308	caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
	309
247abff4	310	if(DebugLevel() > 2)
247abff4	311	{
f2ccb5b8	312	printf("Filter, aod : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
247abff4	313	caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
247abff4	314	caloCluster->GetPosition(posF);
f2ccb5b8	315	printf("Filter, aod : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
247abff4	316	}
247abff4	317
3a58eee6	318	//Matched tracks, just to know if there was any match, the track pointer is useless.
f2ccb5b8	319	//Temporary trick, FIXME
3a58eee6	320	if(bESD){
f2ccb5b8	321	if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990) { //Default value in PHOS 999, in EMCAL 1024, why?
	322	if(DebugLevel() > 2)
	323	printf("* Cluster Track-Matched * dR %f, dZ %f\n",caloCluster->GetEmcCpvDistance(),caloCluster->Chi2());
	324	caloCluster->AddTrackMatched(0x0);
	325	}// fill the array with one entry to signal a possible match
	326	//TArrayI* matchedT = ((AliESDCaloCluster*)cluster)->GetTracksMatched();
	327	//if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
	328	// for (Int_t im = 0; im < matchedT->GetSize(); im++) {
	329	// Int_t iESDtrack = matchedT->At(im);;
	330	// if ((AliVTrack*)InputEvent()->GetTrack(iESDtrack) != 0) {
	331	// caloCluster->AddTrackMatched((AliVTrack*)InputEvent()->GetTrack(iESDtrack));
	332	// }
	333	// }
	334	//}// There is at least a match with a track
3a58eee6	335	}
7a4cf423	336	}
	337	caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
	338	// end of loop on calo clusters
	339
	340	// fill EMCAL cell info
247abff4	341	if ((fCaloFilter==kBoth \|\| fCaloFilter==kEMCAL) && event->GetEMCALCells()) { // protection against missing ESD information
	342	AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
	343	Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
7a4cf423	344
	345	AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
	346	aodEMcells.CreateContainer(nEMcell);
c8fe2783	347	aodEMcells.SetType(AliVCaloCells::kEMCALCell);
247abff4	348	Double_t calibFactor = 1.;
247abff4	349	for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
5ef94e1b	350	Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
	351	fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
	352	fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
	353
247abff4	354	if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn()){
247abff4	355	calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
247abff4	356	}
5ef94e1b	357
	358	if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
	359	aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
	360	//printf("GOOD channel\n");
	361	}
	362	else {
	363	aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
	364	//printf("BAD channel\n");
	365
	366	}
7a4cf423	367	}
	368	aodEMcells.Sort();
	369	}
	370
	371	// fill PHOS cell info
247abff4	372	if ((fCaloFilter==kBoth \|\| fCaloFilter==kPHOS) && event->GetPHOSCells()) { // protection against missing ESD information
	373	AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
	374	Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
7a4cf423	375
	376	AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
	377	aodPHcells.CreateContainer(nPHcell);
c8fe2783	378	aodPHcells.SetType(AliVCaloCells::kPHOSCell);
7a4cf423	379	for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
247abff4	380	aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell));
7a4cf423	381	}
	382	aodPHcells.Sort();
	383	}
	384
	385
	386	return;
	387	}
	388
5ef94e1b	389	//_____________________________________________________
	390	void AliAnalysisTaskCaloFilter::PrintInfo(){
	391
	392	//Print settings
	393
	394	printf("TASK: AnalysisCaloFilter \n");
	395	printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
	396	printf("\t Calorimeter Filtering Option ? %d\n",fCaloFilter);
	397
	398	}
	399
247abff4	400	//_____________________________________________________
	401	//void AliAnalysisTaskCaloFilter::LocalInit()
	402	//{
	403	// // Local Initialization
	404	//
	405	// // Create cuts/param objects and publish to slot
	406	// const Int_t buffersize = 255;
	407	// char onePar[buffersize] ;
	408	// fCuts = new TList();
	409	//
	410	// snprintf(onePar,buffersize, "Calorimeter Filtering Option %d", fCaloFilter) ;
	411	// fCuts->Add(new TObjString(onePar));
	412	// snprintf(onePar,buffersize, "Not only filter but correct? %d cells;", fCorrect) ;
	413	// fCuts->Add(new TObjString(onePar));
	414	//
	415	// // Post Data
	416	// PostData(2, fCuts);
	417	//
	418	//}
7a4cf423	419
7a4cf423	420
	421	//__________________________________________________
	422	void AliAnalysisTaskCaloFilter::Terminate(Option_t /option*/)
	423	{
	424	// Terminate analysis
	425	//
	426	if (fDebug > 1) printf("AnalysisCaloFilter: Terminate() \n");
	427	}
	428