// PHOS or both, creating the corresponing AODCaloClusters
// and AODCaloCells.
// Keep also the AODHeader information and the vertex.
-// Needed for calorimeter calibration.
+// Keep tracks if requested
// Copy of AliAnalysisTaskESDfilter.
// Author: Gustavo Conesa Balbastre (INFN - Frascati)
//////////////////////////////////////////////////////////
//Root
#include "TGeoManager.h"
#include "TFile.h"
-#include "TNtuple.h"
#include "TROOT.h"
#include "TInterpreter.h"
#include "AliAnalysisTaskCaloFilter.h"
ClassImp(AliAnalysisTaskCaloFilter)
-
+
////////////////////////////////////////////////////////////////////////
AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
- AliAnalysisTaskSE("CaloFilterTask"), //fCuts(0x0),
- fCaloFilter(0), fCorrect(kFALSE),
- fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"),
- fEMCALRecoUtils(new AliEMCALRecoUtils),
- fESDtrackCuts(0), fTriggerAnalysis (new AliTriggerAnalysis), fTrackMultEtaCut(0.8),
- fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
- fGeoMatrixSet(kFALSE), fEventNtuple(0),fConfigName(""),fFillAODFile(kTRUE)
+AliAnalysisTaskSE("CaloFilterTask"), //fCuts(0x0),
+fCaloFilter(0), fCorrect(kFALSE),
+fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"),
+fEMCALRecoUtils(new AliEMCALRecoUtils),
+fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
+fGeoMatrixSet(kFALSE),
+fConfigName(""),fFillAODFile(kTRUE), fFillTracks(kFALSE),
+fEnergyCut(0.), fNcellsCut (0), fVzCut(100.)
{
// Default constructor
- fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
+
for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
//for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i] = 0 ;
-
- DefineOutput(1, TNtuple::Class());
}
//__________________________________________________
AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
- AliAnalysisTaskSE(name), //fCuts(0x0),
- fCaloFilter(0), fCorrect(kFALSE),
- fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"),
- fEMCALRecoUtils(new AliEMCALRecoUtils),
- fESDtrackCuts(0), fTriggerAnalysis (new AliTriggerAnalysis), fTrackMultEtaCut(0.8),
- fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
- fGeoMatrixSet(kFALSE), fEventNtuple(0),fConfigName(""),fFillAODFile(kTRUE)
+AliAnalysisTaskSE(name), //fCuts(0x0),
+fCaloFilter(0), fCorrect(kFALSE),
+fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"),
+fEMCALRecoUtils(new AliEMCALRecoUtils),
+fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
+fGeoMatrixSet(kFALSE),
+fConfigName(""),fFillAODFile(kTRUE),fFillTracks(kFALSE),
+fEnergyCut(0.), fNcellsCut (0), fVzCut(100.)
{
// Constructor
- fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
+
for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
//for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i] = 0 ;
-
- DefineOutput(1, TNtuple::Class());
-
+
}
//__________________________________________________
if(fEMCALGeo) delete fEMCALGeo;
if(fEMCALRecoUtils) delete fEMCALRecoUtils;
- if(fESDtrackCuts) delete fESDtrackCuts;
- if(fTriggerAnalysis) delete fTriggerAnalysis;
-
- if(fEventNtuple) delete fEventNtuple;
}
{
//Init analysis with configuration macro if available
- if(gROOT->LoadMacro(fConfigName) >=0){
+ 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;
- fTrackMultEtaCut = filter->fTrackMultEtaCut;
- fESDtrackCuts = filter->fESDtrackCuts;
+ fEnergyCut = filter->fEnergyCut;
+ fNcellsCut = filter->fNcellsCut;
+ fVzCut = filter->fVzCut;
+
for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
}
}
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());
AliVEvent* event = InputEvent();
- if(!event) {
+ if(!event)
+ {
printf("AliAnalysisTaskCaloFilter::UserExec - This event does not contain Input?");
return;
}
-
+
+ // Select the event
+
+ if(!AcceptEventVertex()) return;
+ if(!AcceptEventEMCAL()) return;
+
//Magic line to write events to file
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;
+ 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);
+
+ //
+ FillAODHeader();
+
+ //
+ FillAODVertices();
+
+ //
+ FillAODTracks();
+
+ //
+ CorrectionsInEMCAL();
+
+ //
+ FillAODCaloClusters();
+
+ //
+ FillAODCaloCells();
+
+}
+
+//___________________________________________________
+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
+
+ if(fCaloFilter!=kEMCAL) return kTRUE; // accept all
+
+ Int_t nCluster = InputEvent() -> GetNumberOfCaloClusters();
+ AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
+ Int_t bc = InputEvent() -> GetBunchCrossNumber();
+
+ for(Int_t icalo = 0; icalo < nCluster; icalo++)
+ {
+ AliESDCaloCluster *clus = (AliESDCaloCluster*) (InputEvent()->GetCaloCluster(icalo));
+
+ if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fNcellsCut ) && ( clus->E() > fEnergyCut ) &&
+ fEMCALRecoUtils->IsGoodCluster(clus,fEMCALGeo,caloCell,bc))
+ {
+
+ // printf("Accept event %d, E %2.2f > %2.2f, nCells %d > %d \n",
+ // (Int_t) Entry(),clus->E(), fEnergyCut, clus->GetNCells(), fNcellsCut);
+
+ return kTRUE;
+ }
- }//ESDs
+ }// loop
- 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);
+ return kFALSE;
- //Put bools with event selection parameters in a TNtuple
- //Int_t eventSelection[] = {bPileup,bGoodVertex,bV0AND,trackMult};
- fEventNtuple->Fill(bPileup,bGoodVertex,bV0AND,trackMult);
+}
+
+//________________________________________________
+void AliAnalysisTaskCaloFilter::FillAODCaloCells()
+{
+ // Fill EMCAL/PHOS cell info
- //--------------------------------------------------------------------
- //Set in AOD General Event Parameters, vertex, runnumber, trigger, etc
- //-------------------------------------------------------------------
+ AliVEvent * event = InputEvent();
- // set arrays and pointers
+ // EMCAL
+ 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,
+ 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));
+ }
+
+ aodPHcells.Sort();
+ }
+}
+
+
+//___________________________________________________
+void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
+{
+ // Fill the AOD with caloclusters
+
+ AliVEvent * event = InputEvent();
+
+ // Access to the AOD container of clusters
+
+ TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
+ Int_t jClusters=0;
Float_t posF[3] ;
- Double_t pos[3] ;
- Double_t covVtx[6];
- for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
+
+ Int_t nCaloClus = event->GetNumberOfCaloClusters();
+ 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 ;
+ }
+
+ // 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
+ }
+
+ caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots
+
+}
+
+
+//______________________________________________
+void AliAnalysisTaskCaloFilter::FillAODHeader()
+{
+ // AOD header copy
+
+ AliVEvent* event = InputEvent();
+ AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
+ AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
+
AliAODHeader* header = AODEvent()->GetHeader();
+ // Copy from AODs
+ if(aodevent)
+ {
+ *header = *(aodevent->GetHeader());
+ return;
+ }
+
+ // Copy from ESDs
+
header->SetRunNumber(event->GetRunNumber());
- if(esdevent){
- TTree* tree = fInputHandler->GetTree();
- if (tree) {
- TFile* file = tree->GetCurrentFile();
- if (file) header->SetESDFileName(file->GetName());
- }
+ 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->SetEventType(event->GetEventType());
//Centrality
- if(event->GetCentrality()){
+ if(event->GetCentrality())
+ {
header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
}
else{
//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->SetFiredTriggerClasses(esdevent->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->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());
Float_t diamcov[3];
event->GetDiamondCovXY(diamcov);
header->SetDiamond(diamxy,diamcov);
- if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
- else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
+ header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
- //
- //
- Int_t nVertices = 1 ;/* = prim. vtx*/;
- Int_t nCaloClus = event->GetNumberOfCaloClusters();
+}
+
+//_____________________________________________
+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)
+ {
+ TClonesArray* inTracks = aodevent ->GetTracks();
+ TClonesArray* ouTracks = AODEvent()->GetTracks();
+
+ new (ouTracks) TClonesArray(*inTracks);
+
+ return;
+ }
+
+}
+
+//_______________________________________________
+void AliAnalysisTaskCaloFilter::FillAODVertices()
+{
+ // Copy vertices
+
+ AliVEvent* event = InputEvent();
+ AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
+ AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
+
+ // set arrays and pointers
+ Double_t pos[3] ;
+ Double_t covVtx[6];
+ for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
+
+ // Copy from AODs
+ if(aodevent)
+ {
+ TClonesArray* inVertices = aodevent ->GetVertices();
+ TClonesArray* ouVertices = AODEvent()->GetVertices();
+
+ new (ouVertices) TClonesArray(*inVertices);
+ return;
+ }
- AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
+ // Copy from ESDs
// Access to the AOD container of vertices
- TClonesArray &vertices = *(AODEvent()->GetVertices());
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)
}
AliAODVertex * primary = new(vertices[jVertices++])
- AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
+ 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
- //-------------------------------------------
+}
+
+//__________________________________________________
+void AliAnalysisTaskCaloFilter::CorrectionsInEMCAL()
+{
//If EMCAL, and requested, correct energy, position ...
+
//Need to do this in a separate loop before filling the ESDs because of the track matching recalculations
- if(fCorrect && (fCaloFilter==kEMCAL || fCaloFilter==kBoth) ) {
-
- if(!fGeoMatrixSet){
- if(fLoadEMCALMatrices){
+
+ if(fCorrect && (fCaloFilter==kEMCAL || fCaloFilter==kBoth) )
+ {
+ if(!fGeoMatrixSet)
+ {
+ if(fLoadEMCALMatrices)
+ {
printf("AliAnalysisTaskCaloFilter::UserExec() - Load user defined EMCAL geometry matrices\n");
for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
if(fEMCALMatrix[mod]){
fGeoMatrixSet=kTRUE;
}//SM loop
}//Load matrices
- else if(!gGeoManager){
+ else if(!gGeoManager)
+ {
printf("AliAnalysisTaskCaloFilter::UserExec() - Get geo matrices from data\n");
//Still not implemented in AOD, just a workaround to be able to work at least with ESDs
- if(!strcmp(event->GetName(),"AliAODEvent")) {
+ if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
+ {
if(DebugLevel() > 1)
printf("AliAnalysisTaskCaloFilter Use ideal geometry, values geometry matrix not kept in AODs.\n");
}//AOD
- else {
+ else
+ {
if(DebugLevel() > 1) printf("AliAnalysisTaskCaloFilter Load Misaligned matrices. \n");
- AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
- if(!esd) {
+ AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
+ if(!esd)
+ {
printf("AliAnalysisTaskCaloFilter::UserExec() - This event does not contain ESDs?");
return;
}
- for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
+ for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+ {
//if(DebugLevel() > 1)
esd->GetEMCALMatrix(mod)->Print();
if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
fGeoMatrixSet=kTRUE;
}//ESD
}//Load matrices from Data
-
+
}//first event
//Cluster Loop
- for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
+ Int_t nCaloClus = InputEvent()->GetNumberOfCaloClusters();
+
+ for (Int_t iClust=0; iClust<nCaloClus; ++iClust)
+ {
+ AliVCluster * cluster = InputEvent()->GetCaloCluster(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, distToBad %f\n",iClust,cluster->E(),
}
//Recalculate distance to bad channels, if new list of bad channels provided
- fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, event->GetEMCALCells(), cluster);
-
- if(fEMCALRecoUtils->IsRecalibrationOn()) {
- fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
- fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, event->GetEMCALCells(),cluster);
+ fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, InputEvent()->GetEMCALCells(), cluster);
+
+ if(fEMCALRecoUtils->IsRecalibrationOn())
+ {
+ fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, InputEvent()->GetEMCALCells());
+ fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
fEMCALRecoUtils->RecalculateClusterPID(cluster);
}
-
- fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);
+
+ fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
if(DebugLevel() > 2)
{
cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
}
+
//Recalculate track-matching
- fEMCALRecoUtils->FindMatches(event,0,fEMCALGeo);
+ fEMCALRecoUtils->FindMatches(InputEvent(),0,fEMCALGeo);
} // 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 && 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;
- // }
-
- 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->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(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(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) {
- // 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,
- 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();
- }
-
- // 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),
- eventPHcells.GetTime(iCell),eventPHcells.GetMCLabel(iCell),eventPHcells.GetEFraction(iCell));
- }
- aodPHcells.Sort();
- }
-
- PostData(1, fEventNtuple);
-
- //return;
}
-//____________________________________________________________________________
-Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex(){
+//_______________________________________________________
+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) {
+ // AODs
+ if(!event)
+ {
+ // 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 ;
+
+ return kTRUE;
+ }
+
+ // ESDs
+ if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
+ {
return kTRUE;
}
- if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) {
+ if(event->GetPrimaryVertexTracks()->GetNContributors() < 1)
+ {
// SPD vertex
- if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) {
+ if(event->GetPrimaryVertexSPD()->GetNContributors() > 0)
+ {
//cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
return kTRUE;
}
- if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) {
+ if(event->GetPrimaryVertexSPD()->GetNContributors() < 1)
+ {
// cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
return kFALSE;
}
}
return kFALSE;
-
+
}
-
-//_____________________________________________________
-void AliAnalysisTaskCaloFilter::PrintInfo(){
-
+//_________________________________________
+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);
//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);
-}
-
-//_____________________________________________________
-//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");
+ printf("\t Fill AOD file? %d\n", fFillAODFile);
+ printf("\t Fill Tracks? %d\n" , fFillTracks);
+ printf("Event Selection : EMCAL min E %f, EMCAL NCells %d, Vertex %f\n",fEnergyCut,fNcellsCut,fVzCut);
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff