// Author: Gustavo Conesa Balbastre (INFN - Frascati)
//////////////////////////////////////////////////////////
-#include "AliAnalysisTaskCaloFilter.h"
+//Root
+#include "TGeoManager.h"
+#include "TFile.h"
+#include "TNtuple.h"
+#include "TROOT.h"
+#include "TInterpreter.h"
+
+//STEER
#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"
-#include "TGeoManager.h"
+#include "AliTriggerAnalysis.h"
+
+//EMCAL
+#include "AliEMCALRecoUtils.h"
+#include "AliEMCALGeometry.h"
+
+#include "AliAnalysisTaskCaloFilter.h"
ClassImp(AliAnalysisTaskCaloFilter)
////////////////////////////////////////////////////////////////////////
AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
- AliAnalysisTaskSE(), //fCuts(0x0),
+ AliAnalysisTaskSE("CaloFilterTask"), //fCuts(0x0),
fCaloFilter(0), fCorrect(kFALSE),
fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"),
fEMCALRecoUtils(new AliEMCALRecoUtils),
- fESDtrackCuts(0), fTrackMultEtaCut(0.8),
+ fESDtrackCuts(0), fTriggerAnalysis (new AliTriggerAnalysis), fTrackMultEtaCut(0.8),
fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
- fGeoMatrixSet(kFALSE)
+ fGeoMatrixSet(kFALSE), fEventNtuple(0),fConfigName(""),fFillAODFile(kTRUE)
{
// Default constructor
fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
//for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i] = 0 ;
+ DefineOutput(1, TNtuple::Class());
+
}
//__________________________________________________
fCaloFilter(0), fCorrect(kFALSE),
fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"),
fEMCALRecoUtils(new AliEMCALRecoUtils),
- fESDtrackCuts(0), fTrackMultEtaCut(0.8),
+ fESDtrackCuts(0), fTriggerAnalysis (new AliTriggerAnalysis), fTrackMultEtaCut(0.8),
fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
- fGeoMatrixSet(kFALSE)
+ fGeoMatrixSet(kFALSE), fEventNtuple(0),fConfigName(""),fFillAODFile(kTRUE)
{
// Constructor
fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
//for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i] = 0 ;
+ DefineOutput(1, TNtuple::Class());
+
}
//__________________________________________________
{
//Destructor.
- if(fEMCALGeo) delete fEMCALGeo;
- if(fEMCALRecoUtils) delete fEMCALRecoUtils;
- if(fESDtrackCuts) delete fESDtrackCuts;
+ if(fEMCALGeo) delete fEMCALGeo;
+ if(fEMCALRecoUtils) delete fEMCALRecoUtils;
+ if(fESDtrackCuts) delete fESDtrackCuts;
+ if(fTriggerAnalysis) delete fTriggerAnalysis;
+ if(fEventNtuple) delete fEventNtuple;
+
}
+//-------------------------------------------------------------------
+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;
+ fEMCALRecoUtils = filter->fEMCALRecoUtils;
+ fConfigName = filter->fConfigName;
+ fCaloFilter = filter->fCaloFilter;
+ fCorrect = filter->fCorrect;
+ fTrackMultEtaCut = filter->fTrackMultEtaCut;
+ fESDtrackCuts = filter->fESDtrackCuts;
+ for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
+ }
+}
+
//__________________________________________________
void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
{
fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
+ OpenFile(1);
+
+ fEventNtuple = new TNtuple("EventSelection","EventSelection", "bPileup:bGoodVertex:bV0AND:trackMult");
+
+ PostData(1, fEventNtuple);
+
}
//__________________________________________________
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());
- // 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?");
+ printf("AliAnalysisTaskCaloFilter::UserExec - 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){
+ AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
+
+ // cast event, depending on input we will have to use one or the other type of event
+ AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
+ AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
+
+ //-------------------------------------------
+ //Event selection parameters
+ //-------------------------------------------
+ //Is it a pileup event?
+ Bool_t bPileup = event->IsPileupFromSPD(3, 0.8, 3., 2., 5.); //Default values, if not it does not compile
+ //Bool_t bPileup = event->IsPileupFromSPD();
+ //if(bPileup) return kFALSE;
+ Int_t trackMult = 0;
+ Bool_t bV0AND = kFALSE;
+ Bool_t bGoodVertex = kFALSE;
+ if(esdevent){
+ //Get track multiplicity
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.8
if(TMath::Abs(track->Eta())< fTrackMultEtaCut) trackMult++;
- }
- }
+ }
+ //V0AND?
+ bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esdevent, AliTriggerAnalysis::kV0AND);
+ //if(!bV0AND) return kFALSE;
+ //Well reconstructed vertex
+ bGoodVertex = CheckForPrimaryVertex();
+ //if(!bGoodVertex) return kFALSE;
+
+ }//ESDs
- // set arrays and pointers
- Float_t posF[3];
- Double_t pos[3];
+ if(fDebug > 0)
+ printf("AliAnalysisTaskCaloFilter::UserExec() - PileUp %d, Good Vertex %d, v0AND %d, Track Mult in |eta| < %2.1f = %d\n",
+ bPileup,bGoodVertex,bV0AND, fTrackMultEtaCut, trackMult);
- Double_t covVtx[6];
+ //Put bools with event selection parameters in a TNtuple
+ //Int_t eventSelection[] = {bPileup,bGoodVertex,bV0AND,trackMult};
+ fEventNtuple->Fill(bPileup,bGoodVertex,bV0AND,trackMult);
+
+ //--------------------------------------------------------------------
+ //Set in AOD General Event Parameters, vertex, runnumber, trigger, etc
+ //-------------------------------------------------------------------
+ // 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());
+
+ if(esdevent){
+ TTree* tree = fInputHandler->GetTree();
+ if (tree) {
+ TFile* file = tree->GetCurrentFile();
+ if (file) header->SetESDFileName(file->GetName());
+ }
+ }
+ else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
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());
+ //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
+ if (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
+ else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
header->SetTriggerMask(event->GetTriggerMask());
header->SetTriggerCluster(event->GetTriggerCluster());
+ if(esdevent){
+ header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
+ header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
+ header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
+ }
+ else if (aodevent){
+ header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
+ header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
+ header->SetL2TriggerInputs(aodevent->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);
+ Float_t diamcov[3];
+ event->GetDiamondCovXY(diamcov);
header->SetDiamond(diamxy,diamcov);
- if(bESD){
- header->SetDiamondZ(((AliESDEvent*)event)->GetDiamondZ(),((AliESDEvent*)event)->GetSigma2DiamondZ());
- }
+ if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
+ else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
+
//
//
Int_t nVertices = 1 ;/* = prim. vtx*/;
// 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();
+ if (esdevent){
+ esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
+ chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
}
- else if (bAOD){
- ((AliAODEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
- chi = ((AliAODEvent*)event)->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
+ else if (aodevent){
+ aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
+ chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
}
AliAODVertex * primary = new(vertices[jVertices++])
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) ) {
if(!fGeoMatrixSet){
if(fLoadEMCALMatrices){
- printf("AliAnalysisTaskCaloFilter::UserExec() - Load user defined geometry matrices\n");
+ printf("AliAnalysisTaskCaloFilter::UserExec() - Load user defined EMCAL geometry matrices\n");
for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
if(fEMCALMatrix[mod]){
if(DebugLevel() > 1)
fGeoMatrixSet=kTRUE;
}//ESD
}//Load matrices from Data
+
+ //Recover time dependent corrections, put then in recalibration histograms. Do it once
+ fEMCALRecoUtils->SetRunDependentCorrections(InputEvent()->GetRunNumber());
+
}//first event
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, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
}
+ cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
+
}
//Recalculate track-matching
- fEMCALRecoUtils->FindMatches(event);
+ fEMCALRecoUtils->FindMatches(event,0,fEMCALGeo);
} // corrections in EMCAL
- //Now loop on clusters to fill AODs
+ //-------------------------------------------
+ // Now loop on clusters to fill AODs
+ //-------------------------------------------
for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
AliVCluster * cluster = event->GetCaloCluster(iClust);
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(cluster->IsEMCAL() && fCorrect && esdevent){
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);
+ 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
AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
AliAODCaloCluster(id,
- 0,
- 0x0,
+ cluster->GetNLabels(),
+ cluster->GetLabels(),
energy,
posF,
NULL,
//Matched tracks, just to know if there was any match, the track pointer is useless.
//Temporary trick, FIXME
- if(bESD){
+ if(esdevent){
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);
+ caloCluster->AddTrackMatched(new AliAODTrack);
}// fill the array with one entry to signal a possible match
//TArrayI* matchedT = ((AliESDCaloCluster*)cluster)->GetTracksMatched();
//if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
else {
aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
//printf("BAD channel\n");
-
}
}
aodEMcells.Sort();
aodPHcells.Sort();
}
+ PostData(1, fEventNtuple);
- return;
+ //return;
}
+//____________________________________________________________________________
+Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex(){
+ //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
+ //It only works for ESDs
+
+ AliESDEvent * event = dynamic_cast<AliESDEvent*> (InputEvent());
+ if(!event) return kFALSE;
+
+ if(event->GetPrimaryVertexTracks()->GetNContributors() > 0) {
+ return kTRUE;
+ }
+
+ if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) {
+ // SPD vertex
+ if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) {
+ //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+ return kTRUE;
+
+ }
+ if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) {
+ // cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+ return kFALSE;
+ }
+ }
+ return kFALSE;
+
+}
+
+
//_____________________________________________________
void AliAnalysisTaskCaloFilter::PrintInfo(){
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);
-
}
//_____________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff