#include "AliMCEvent.h"
#include "AliAODMCHeader.h"
#include "AliGenPythiaEventHeader.h"
+#include "AliVEvent.h"
#include "AliAODEvent.h"
+#include "AliVTrack.h"
+#include "AliVParticle.h"
+#include "AliMixedEvent.h"
ClassImp(AliCaloTrackReader)
fAODPHOSNormalInputEntries(0), fTrackStatus(0),
fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
fCleanOutputStdAOD(kFALSE), fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
- fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0)
+ fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0),
+ fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL), fWriteOutputStdAOD(kFALSE)
{
//Ctor
fPHOSCells->Clear() ;
delete fPHOSCells ;
}
+ if (fMixedEvent) {
+ for (Int_t i = 0; i < fNMixedEvent; i++) {
+ delete [] fVertex[i] ;
+ }
+ }
+ delete [] fVertex ;
// Pointers not owned, done by the analysis frame
// if(fInputEvent) delete fInputEvent ;
//____________________________________________________________________________
TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const {
- //Return list of particles in AOD. Do it for the corresponding input event.
- if(fDataType == kAOD){
- //Normal input AOD
- if(input == 0) return (TClonesArray*)((AliAODEvent*)fInputEvent)->FindListObject("mcparticles");
- //Second input AOD
- else if(input == 1 && fSecondInputAODEvent) return (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");
- else {
- printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input);
- return 0x0;
- }
- }
- else {
- printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
- return 0x0;
- }
+ //Return list of particles in AOD. Do it for the corresponding input event.
+ TClonesArray * rv = NULL ;
+ // if(fDataType == kAOD){
+ //Normal input AOD
+ AliAODEvent * evt = NULL ;
+ if (evt) {
+ if(input == 0){
+ evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
+ rv = (TClonesArray*)evt->FindListObject("mcparticles");
+ } else if(input == 1 && fSecondInputAODEvent){ //Second input AOD
+ rv = (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");
+ } else {
+ printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input);
+ }
+ } else {
+ printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
+ }
+ return rv ;
}
//____________________________________________________________________________
if(fInputEvent->GetHeader())
eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
if( fFiredTriggerClassName !="" && !fAnaLED){
- if(eventType!=7)return kFALSE; //Only physics event, do not use for simulated events!!!
+ if(eventType!=7)
+ return kFALSE; //Only physics event, do not use for simulated events!!!
if(fDebug > 0)
printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
}
- if(fFillEMCALCells) FillInputEMCALCells();
- if(fFillPHOSCells) FillInputPHOSCells();
+
+ for (Int_t iev = 0; iev < fNMixedEvent; iev++) {
+ if (!fMixedEvent) {
+ GetVertex() ;
+ }
+ else
+ fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
+ }
+
+ if(fFillEMCALCells)
+ FillInputEMCALCells();
+ if(fFillPHOSCells)
+ FillInputPHOSCells();
- if(fFillCTS) FillInputCTS();
- if(fFillEMCAL) FillInputEMCAL();
- if(fFillPHOS) FillInputPHOS();
+ if(fFillCTS)
+ FillInputCTS();
+ if(fFillEMCAL)
+ FillInputEMCAL();
+ if(fFillPHOS)
+ FillInputPHOS();
return kTRUE ;
}
}
+
+//____________________________________________________________________________
+void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
+{
+ fInputEvent = input;
+ fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
+ if (fMixedEvent) {
+ fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
+ }
+ fVertex = new Double_t*[fNMixedEvent] ;
+ for (Int_t i = 0; i < fNMixedEvent; i++) {
+ fVertex[i] = new Double_t[3] ;
+ fVertex[i][0] = 0.0 ;
+ fVertex[i][1] = 0.0 ;
+ fVertex[i][2] = 0.0 ;
+ }
+}
+
+//____________________________________________________________________________
+Double_t * AliCaloTrackReader::GetVertex() {
+ //Return vertex position
+ if (fMixedEvent)
+ return NULL ;
+ if(fInputEvent->GetPrimaryVertex())
+ fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
+ else {
+ printf("AliCaloTrackReader::GetVertex() - No vertex available, InputEvent()->GetPrimaryVertex() = 0, return (0,0,0)\n");
+ fVertex[0][0]=0.0; fVertex[0][1]=0.0; fVertex[0][2]=0.0;
+ }
+ return fVertex[0] ;
+}
+
+ //____________________________________________________________________________
+void AliCaloTrackReader::FillInputCTS() {
+ //Return array with Central Tracking System (CTS) tracks
+ if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputCTS()\n");
+ Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
+ Int_t naod = 0;
+ Double_t p[3];
+ for (Int_t itrack = 0; itrack < nTracks; itrack++) {////////////// track loop
+ AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
+
+ //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
+ if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus))
+ continue ;
+
+ track->GetPxPyPz(p) ;
+ TLorentzVector momentum(p[0],p[1],p[2],0);
+
+ if(fCTSPtMin < momentum.Pt()){
+
+ if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS"))
+ continue;
+
+ if(fDebug > 2 && momentum.Pt() > 0.1)
+ printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
+ momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
+
+ if (fMixedEvent) {
+ track->SetID(itrack);
+ }
+ if(fWriteOutputStdAOD){
+ AliAODTrack * aodTrack = dynamic_cast<AliAODTrack*> (track) ;
+ if (aodTrack) {
+ AliAODTrack* newtrack = new((*(fOutputEvent->GetTracks()))[naod++]) AliAODTrack(*aodTrack);
+ fAODCTS->Add(newtrack); //Use AOD stored in output for references.
+ } else {
+ AliError("Can only write AOD tracks to output AOD !") ;
+ }
+ }
+ else {
+ fAODCTS->Add(track);
+ }
+ }//Pt and Fiducial cut passed.
+ }// track loop
+
+ fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
+ if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS() - aod entries %d\n", fAODCTSNormalInputEntries);
+
+ // //If second input event available, add the clusters.
+ // if(fSecondInputAODTree && fSecondInputAODEvent){
+ // nTracks = fSecondInputAODEvent->GetNumberOfTracks() ;
+ // if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS() - Add second input tracks, entries %d\n", nTracks) ;
+ // for (Int_t itrack = 0; itrack < nTracks; itrack++) {////////////// track loop
+ // AliAODTrack * track = ((AliAODEvent*)fSecondInputAODEvent)->GetTrack(itrack) ; // retrieve track from esd
+ //
+ // //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
+ // if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
+ //
+ // track->GetPxPyPz(p) ;
+ // TLorentzVector momentum(p[0],p[1],p[2],0);
+ //
+ // if(fCTSPtMin < momentum.Pt()){
+ //
+ // if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
+ //
+ // if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
+ // momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
+ //
+ // if(fWriteOutputStdAOD){
+ // AliAODTrack* newtrack = new((*(fOutputEvent->GetTracks()))[naod++]) AliAODTrack(*track);
+ // fAODCTS->Add(newtrack); //Use AOD stored in output for references.
+ // }
+ // else fAODCTS->Add(track);
+ //
+ // }//Pt and Fiducial cut passed.
+ // }// track loop
+ //
+ // if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS() - aod normal entries %d, after second input %d\n", fAODCTSNormalInputEntries, fAODCTS->GetEntriesFast());
+ // } //second input loop
+ //
+}
+
+ //____________________________________________________________________________
+void AliCaloTrackReader::FillInputEMCAL() {
+ //Return array with EMCAL clusters in aod format
+ if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputEMCAL()\n");
+
+ Int_t naod = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
+ //Loop to select clusters in fiducial cut and fill container with aodClusters
+ Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
+ for (Int_t iclus = 0; iclus < nclusters; iclus++) {
+ AliVCluster * clus = 0;
+ if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
+ if (clus->IsEMCAL()){
+
+ //Check if the cluster contains any bad channel and if close to calorimeter borders
+ Int_t vindex = 0 ;
+ if (fMixedEvent)
+ vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
+
+ if(GetCaloUtils()->ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells()))
+ continue;
+ if(!GetCaloUtils()->CheckCellFiducialRegion(clus, (AliVCaloCells*)fInputEvent->GetEMCALCells(), fInputEvent, vindex))
+ continue;
+
+ TLorentzVector momentum ;
+
+ clus->GetMomentum(momentum, fVertex[vindex]);
+
+ if(fEMCALPtMin < momentum.Pt()){
+
+ if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL"))
+ continue;
+
+ if(fDebug > 2 && momentum.E() > 0.1)
+ printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
+ momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
+
+ //Recalibrate the cluster energy
+ if(GetCaloUtils()->IsRecalibrationOn()) {
+ Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetEMCALCells());
+ clus->SetE(energy);
+ }
+
+ if (fMixedEvent) {
+ clus->SetID(iclus) ;
+ }
+ if(fWriteOutputStdAOD){
+ AliAODCaloCluster * aodClus = dynamic_cast<AliAODCaloCluster*> (clus) ;
+ if (aodClus) {
+ AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*aodClus);
+ fAODEMCAL->Add(newclus);
+ } else {
+ AliError("Can only write AOD clusters to output AOD !") ;
+ }
+ } else {
+ fAODEMCAL->Add(clus);
+ }
+ }//Pt and Fiducial cut passed.
+ }//EMCAL cluster
+ }// cluster exists
+ }// cluster loop
+ fAODEMCALNormalInputEntries = fAODEMCAL->GetEntriesFast();
+ if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod entries %d\n", fAODEMCALNormalInputEntries);
+
+ //If second input event available, add the clusters.
+ // if(fSecondInputAODTree && fSecondInputAODEvent){
+ // GetSecondInputAODVertex(v);
+ // nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNumberOfCaloClusters();
+ // if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Add second input clusters, entries %d\n", nclusters) ;
+ // for (Int_t iclus = 0; iclus < nclusters; iclus++) {
+ // AliAODCaloCluster * clus = 0;
+ // if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
+ // if (clus->IsEMCAL()){
+ // TLorentzVector momentum ;
+ // clus->GetMomentum(momentum, v);
+ //
+ // if(fEMCALPtMin < momentum.Pt()){
+ //
+ // if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) continue;
+ //
+ // if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
+ // momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
+ // if(fWriteOutputStdAOD){
+ // AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
+ // fAODEMCAL->Add(newclus);
+ // }
+ // else fAODEMCAL->Add(clus);
+ // }//Pt and Fiducial cut passed.
+ // }//EMCAL cluster
+ // }// cluster exists
+ // }// cluster loop
+ //
+ // if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fAODEMCALNormalInputEntries, fAODEMCAL->GetEntriesFast());
+ //
+ // } //second input loop
+}
+
+ //____________________________________________________________________________
+void AliCaloTrackReader::FillInputPHOS() {
+ //Return array with PHOS clusters in aod format
+ if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputPHOS()\n");
+
+ //Get vertex for momentum calculation
+
+ Int_t naod = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
+ //Loop to select clusters in fiducial cut and fill container with aodClusters
+ Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
+ for (Int_t iclus = 0; iclus < nclusters; iclus++) {
+ AliVCluster * clus = 0;
+ if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
+ if (clus->IsPHOS()){
+ //Check if the cluster contains any bad channel and if close to calorimeter borders
+ Int_t vindex = 0 ;
+ if (fMixedEvent)
+ vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
+ if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells()))
+ continue;
+ if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex))
+ continue;
+
+ TLorentzVector momentum ;
+
+ clus->GetMomentum(momentum, fVertex[vindex]);
+
+ if(fPHOSPtMin < momentum.Pt()){
+
+ if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS"))
+ continue;
+
+ if(fDebug > 2 && momentum.E() > 0.1)
+ printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
+ momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
+
+ //Recalibrate the cluster energy
+ if(GetCaloUtils()->IsRecalibrationOn()) {
+ Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
+ clus->SetE(energy);
+ }
+
+ if (fMixedEvent) {
+ clus->SetID(iclus) ;
+ }
+
+ if(fWriteOutputStdAOD){
+ AliAODCaloCluster * aodClus = dynamic_cast<AliAODCaloCluster*> (clus) ;
+ if (aodClus) {
+ AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*aodClus);
+ fAODPHOS->Add(newclus);
+ } else {
+ AliError("Can only write AOD clusters to output AOD !") ;
+ }
+ } else {
+ fAODPHOS->Add(clus);
+ }
+ }//Pt and Fiducial cut passed.
+ }//PHOS cluster
+ }//cluster exists
+ }//esd cluster loop
+
+ fAODPHOSNormalInputEntries = fAODPHOS->GetEntriesFast() ;
+ if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS() - aod entries %d\n", fAODPHOSNormalInputEntries);
+
+ //If second input event available, add the clusters.
+ // if(fSecondInputAODTree && fSecondInputAODEvent){
+ // GetSecondInputAODVertex(v);
+ // nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNumberOfCaloClusters();
+ // if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS() - Add second input clusters, entries %d\n", nclusters);
+ // for (Int_t iclus = 0; iclus < nclusters; iclus++) {
+ // AliAODCaloCluster * clus = 0;
+ // if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
+ // if (clus->IsPHOS()){
+ // TLorentzVector momentum ;
+ // clus->GetMomentum(momentum, v);
+ //
+ // if(fPHOSPtMin < momentum.Pt()){
+ //
+ // if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
+ //
+ // if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
+ // momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
+ // if(fWriteOutputStdAOD){
+ // AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
+ // fAODPHOS->Add(newclus);
+ // }
+ // else fAODPHOS->Add(clus);
+ // }//Pt and Fiducial cut passed.
+ // }//PHOS cluster
+ // }// cluster exists
+ // }// cluster loop
+ // if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS() - aod normal entries %d, after second input %d\n", fAODPHOSNormalInputEntries, fAODPHOS->GetEntriesFast());
+ // } //second input loop
+
+}
+
+ //____________________________________________________________________________
+void AliCaloTrackReader::FillInputEMCALCells() {
+ //Return array with EMCAL cells in aod format
+
+ fEMCALCells = fInputEvent->GetEMCALCells();
+
+}
+
+ //____________________________________________________________________________
+void AliCaloTrackReader::FillInputPHOSCells() {
+ //Return array with PHOS cells in aod format
+
+ fPHOSCells = fInputEvent->GetPHOSCells();
+
+}
+