#include "AliMCEvent.h"
#include "AliStack.h"
#include "TParticle.h"
+#include "AliAODMCParticle.h"
#include "AliESDtrackCuts.h"
#include "AliV0vertexer.h"
#include "AliESDCaloCluster.h"
#include "AliESDCaloCells.h"
+#include "AliAODEvent.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALRecoUtils.h"
+#include "AliOADBContainer.h"
+#include "AliVEvent.h"
+#include "AliVVertex.h"
#include "AliVCluster.h"
+#include "AliVCaloCells.h"
#include "AliAnalysisTaskEMCALClusterizeFast.h"
#include "TLorentzVector.h"
fPrTrCuts(0),
fSelTracks(0),
fSelPrimTracks(0),
+ fTracks(0),
fPhotConvArray(0),
fMyClusts(0),
fMyAltClusts(0),
fMyTracks(0),
fMyMcParts(0),
fHeader(0x0),
+ fOADBContainer(0),
fCaloClusters(0),
fCaloClustersNew(0),
- fEMCalCells(0),
+ fAODMCParticles(0),
+ fVCells(0),
fGeom(0x0),
fTimeResTOF(0),
fMipResponseTPC(0),
fIsTrain(0),
fIsMC(0),
fDebug(0),
+ fRedoV0(1),
fIsGrid(0),
fClusThresh(2.0),
fClusterizer(0),
fCaloClustersName("EmcalClusterv2"),
fESD(0),
+ fAOD(0),
+ fVev(0),
fMCEvent(0),
fStack(0x0),
fOutputList(0),
fTree(0),
+ fMyMcIndex(0),
fNV0sBefAndAftRerun(0),
fConversionVtxXY(0),
fInvMassV0(0),
fDedxPAll(0)
{
// Default constructor.
+ for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
}
//________________________________________________________________________
fPrTrCuts(0),
fSelTracks(0),
fSelPrimTracks(0),
+ fTracks(0),
fPhotConvArray(0),
fMyClusts(0),
fMyAltClusts(0),
fMyTracks(0),
fMyMcParts(0),
fHeader(0),
+ fOADBContainer(0),
fCaloClusters(0),
fCaloClustersNew(0),
- fEMCalCells(0),
+ fAODMCParticles(0),
+ fVCells(0),
fGeom(0x0),
fTimeResTOF(0),
fMipResponseTPC(0),
fIsTrain(0),
fIsMC(0),
fDebug(0),
+ fRedoV0(1),
fIsGrid(0),
fClusThresh(2.),
fClusterizer(0),
fCaloClustersName("EmcalClusterv2"),
fESD(0),
+ fAOD(0),
+ fVev(0),
fMCEvent(0),
fStack(0x0),
fOutputList(0),
fTree(0),
+ fMyMcIndex(0),
fNV0sBefAndAftRerun(0),
fConversionVtxXY(0),
fInvMassV0(0),
void AliAnalysisTaskEMCALPhoton::UserCreateOutputObjects()
{
// Create histograms, called once.
-
+ if(this->fDebug)
+ printf("::UserCreateOutputObjects() starting\n");
+
fSelTracks = new TObjArray();
fSelPrimTracks = new TObjArray();
fMyMcParts = new TClonesArray("AliPhotonMcPartObj");
}
- fCaloClusters = new TRefArray();
+ fCaloClusters = new TClonesArray();
fOutputList = new TList();
fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
if( !fTree){
- TFile *f = OpenFile(1);
+ TFile *f = OpenFile(2);
TDirectory::TContext context(f);
if (f) {
f->SetCompressionLevel(2);
fTree->Branch("mcparts", &fMyMcParts, 8*16*1024, 99);
}
}
- if(fIsGrid)fOutputList->Add(fTree);
+ //if(fIsGrid)fOutputList->Add(fTree);
fGeom = AliEMCALGeometry::GetInstance(fGeoName);
+ fOADBContainer = new AliOADBContainer("AliEMCALgeo");
+ fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
fNV0sBefAndAftRerun = new TH2F("hNV0sBefAndAftRerun","check if the number of v0s change with rerun;old v0 n;new v0 n",50,0.5,50.5,50,0.5,50.5);
PostData(1, fOutputList);
- // PostData(2, fTree);
+ PostData(2, fTree);
+
+ if(this->fDebug)
+ printf("::UserCreateOutputObjects() DONE!\n");
+
}
//________________________________________________________________________
{
// User exec, called once per event.
- Bool_t isSelected = 0;
+ Bool_t isSelected = kTRUE;
if(fPeriod.Contains("11")){
if(fPeriod.Contains("11a"))
isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
if(fPeriod.Contains("11c") ||fPeriod.Contains("11d") )
isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
if(fPeriod.Contains("11h") )
- isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA);
+ isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);//kEMCEGA);
}
if(fIsMC){
isSelected = kTRUE;
+ if(this->fDebug)
+ printf("+++Message+++: MC input files.\n");
+ }
+ if(!isSelected){
+ if(this->fDebug)
+ printf("+++Message+++: Event didn't pass the selection\n");
+ return;
+ }
+ if(this->fDebug){
+ printf("::UserExec(): event accepted\n");
+ if(fIsMC)
+ printf("\t in MC mode\n");
}
+ TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
+ TFile *inpfile = (TFile*)tree->GetCurrentFile();
- // Post output data.
- fESD = dynamic_cast<AliESDEvent*>(InputEvent());
- if (!fESD) {
- printf("ERROR: fESD not available, returning...\n");
+ fVev = (AliVEvent*)InputEvent();
+ if (!fVev) {
+ printf("ERROR: event not available\n");
return;
}
+ Int_t runnumber = InputEvent()->GetRunNumber() ;
+ if(this->fDebug)
+ printf("run number = %d\n",runnumber);
+ fESD = dynamic_cast<AliESDEvent*>(fVev);
+ if(!fESD){
+ fAOD = dynamic_cast<AliAODEvent*>(fVev);
+ if(!fAOD){
+ printf("ERROR: Invalid type of event!!!\n");
+ return;
+ }
+ else if(this->fDebug)
+ printf("AOD EVENT!\n");
+ }
- AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
- if(!pv)
+ AliVVertex *pv = (AliVVertex*)fVev->GetPrimaryVertex();
+ Bool_t pvStatus = kTRUE;
+ if(fESD){
+ AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex();
+ pvStatus = esdv->GetStatus();
+ }
+
+ if(!pv) {
+ printf("Error: no primary vertex found!\n");
return;
+ }
+ if(!pvStatus && this->fDebug)
+ printf("bad pv status\n");
if(TMath::Abs(pv->GetZ())>15)
return;
+ if(this->fDebug)
+ printf("+++Message+++: event passed the vertex cut\n");
- TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
- TFile *inpfile = (TFile*)tree->GetCurrentFile();
+ if(fESD)
+ fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
+ if(fAOD)
+ fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
+
+ if(!fTracks){
+ AliError("Track array in event is NULL!");
+ if(this->fDebug)
+ printf("returning due to not finding tracks!\n");
+ return;
+ }
+ const Int_t Ntracks = fTracks->GetEntriesFast();
// Track loop to fill a pT spectrum
- for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
- AliESDtrack* track = fESD->GetTrack(iTracks);
+ for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
+ AliVTrack *track = (AliVTrack*)fTracks->At(iTracks);
if (!track)
continue;
-
-
- if (fTrCuts && fTrCuts->IsSelected(track)){
- fSelTracks->Add(track);
- fDedxPAll->Fill(track->P(), track->GetTPCsignal());
+ AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
+ AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
+ if (esdTrack){
+ if (fTrCuts && fTrCuts->IsSelected(track)){
+ fSelTracks->Add(track);
+ fDedxPAll->Fill(track->P(), track->GetTPCsignal());
+ }
+ if (fPrTrCuts && fPrTrCuts->IsSelected(track))
+ fSelPrimTracks->Add(track);
}
- if (fPrTrCuts && fPrTrCuts->IsSelected(track))
+ else if(aodTrack)
fSelPrimTracks->Add(track);
- } //track loop
+ }//track loop
+
+
fHeader->fInputFileName = inpfile->GetName();
- fHeader->fTrClassMask = fESD->GetHeader()->GetTriggerMask();
- fHeader->fTrCluster = fESD->GetHeader()->GetTriggerCluster();
+ fHeader->fRunNumber = runnumber;
+ fHeader->fTrClassMask = fVev->GetHeader()->GetTriggerMask();
+ fHeader->fTrCluster = fVev->GetHeader()->GetTriggerCluster();
AliCentrality *cent = InputEvent()->GetCentrality();
fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
- if(!fIsTrain){
- for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
- if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
- break;
+
+ TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
+ for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
+ if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
+ break;
+ /*if(fESD)
fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
- }
+ else*/
+ // if(event->GetEMCALMatrix(mod))
+ fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
+ fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
+ }
+ Int_t trackMult = 0;
+ if(fESD){
+ AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
+ trackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
+ if(this->fDebug)
+ printf("ESD Track mult= %d\n",trackMult);
+ }
+ else if(fAOD){
+ trackMult = Ntracks;
+ if(this->fDebug)
+ printf("AOD Track mult= %d\n",trackMult);
+ }
+ if(fESD){
+ TList *l = fESD->GetList();
+ fCaloClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
+ fVCells = fESD->GetEMCALCells();
+ fHeader->fNCells = fVCells->GetNumberOfCells();
+ if(this->fDebug)
+ printf("ESD cluster mult= %d\n",fCaloClusters->GetEntriesFast());
}
- fESD->GetEMCALClusters(fCaloClusters);
- fHeader->fNClus = fCaloClusters->GetEntries();
- fEMCalCells = fESD->GetEMCALCells();
- fHeader->fNCells = fEMCalCells->GetNumberOfCells();
- AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
- fHeader->fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
+ else if(fAOD){
+ fCaloClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
+ fVCells = fAOD->GetEMCALCells();
+ fHeader->fNCells = fVCells->GetNumberOfCells();
+ if(this->fDebug)
+ printf("AOD cluster mult= %d\n",fCaloClusters->GetEntriesFast());
+ }
+
+
+ fHeader->fNClus = fCaloClusters->GetEntriesFast();
+ fHeader->fTrackMult = trackMult;
fMCEvent = MCEvent();
- if(fMCEvent)
+ if(fMCEvent){
fStack = (AliStack*)fMCEvent->Stack();
+ if(this->fStack)
+ fHeader->fNMcParts = this->fStack->GetNtrack();
+ else{
+ printf("Stack not found\n");
+ fHeader->fNMcParts = 0;
+ fAODMCParticles = (TClonesArray*)fVev->FindListObject("mcparticles");
+ }
+ if(fAODMCParticles)
+ fHeader->fNMcParts = fAODMCParticles->GetEntriesFast();
+ }
+ else{
+ if(fIsMC){
+ printf("ERROR: MC Event not available, returning...\n");
+ return;
+ }
+ }
FindConversions();
+ if(this->fDebug)
+ printf("FindConversions done\n");
FillMyCells();
+ if(this->fDebug)
+ printf("FillMyCells done\n");
FillMyClusters();
+ if(this->fDebug)
+ printf("FillMyClusters done\n");
if(fCaloClustersNew)
FillMyAltClusters();
FillIsoTracks();
if(fIsMC)
GetMcParts();
+
+ if(this->fDebug && fIsMC)
+ printf("fMyMcParts nentries=%d",fMyMcParts->GetEntries());
fTree->Fill();
fSelTracks->Clear();
fMyTracks->Clear();
if(fMyMcParts)
fMyMcParts->Clear();
+ fMyMcIndex = 0;
fCaloClusters->Clear();
if(fCaloClustersNew)
fCaloClustersNew->Clear();
+ if(fVCells)
+ fVCells = 0;
+ // Post output data.
PostData(1, fOutputList);
- //PostData(2, fTree);
+ PostData(2, fTree);
}
//________________________________________________________________________
-void AliAnalysisTaskEMCALPhoton::FindConversions()
+void AliAnalysisTaskEMCALPhoton::FindConversions() //WARNING: not ready to use with AODs
{
// Find conversion.
-
+ if(!fESD)//not working with AODs yet
+ return;
+ if(this->fDebug)
+ printf("::FindConversions() starting\n");
if(!fTrCuts)
return;
Int_t iconv = 0;
- Int_t nV0Orig = fESD->GetNumberOfV0s();
- Int_t ntracks = fESD->GetNumberOfTracks();
- fESD->ResetV0s();
- AliV0vertexer lV0vtxer;
- lV0vtxer.Tracks2V0vertices(fESD);
- Int_t nV0s = fESD->GetNumberOfV0s();
+ Int_t nV0Orig = 0;
+ if(fESD)
+ nV0Orig = fESD->GetNumberOfV0s();
+ if(fAOD)
+ nV0Orig = fAOD->GetNumberOfV0s();
+ Int_t nV0s = nV0Orig;
+ Int_t ntracks = fSelTracks->GetEntriesFast();
+ if(fRedoV0 && !fAOD){
+ fESD->ResetV0s();
+ AliV0vertexer lV0vtxer;
+ lV0vtxer.Tracks2V0vertices(fESD);
+ nV0s = fESD->GetNumberOfV0s();
+ }
fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
for(Int_t iv0=0; iv0<nV0s; iv0++){
AliESDv0 * v0 = fESD->GetV0(iv0);
continue;
if(ineg<0 || ineg> ntracks)
continue;
- AliESDtrack *pos = static_cast<AliESDtrack*>(fESD->GetTrack(ipos));
- AliESDtrack *neg = static_cast<AliESDtrack*>(fESD->GetTrack(ineg));
+ AliESDtrack *pos = static_cast<AliESDtrack*>(fSelTracks->At(ipos));
+ AliESDtrack *neg = static_cast<AliESDtrack*>(fSelTracks->At(ineg));
if(!pos || !neg)
continue;
- if (!fTrCuts->IsSelected(pos) || !fTrCuts->IsSelected(neg))
- continue;
/*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
continue;*/
const AliExternalTrackParam * paramPos = v0->GetParamP() ;
myconv->fPosDedx = pos->GetTPCsignal();
myconv->fPosMcLabel = pos->GetLabel();
}
+ if(this->fDebug)
+ printf("::FindConversions() returning...\n\n");
return;
}
void AliAnalysisTaskEMCALPhoton::FillMyCells()
{
// Fill cells.
+ if(this->fDebug)
+ printf("::FillMyCells() starting\n");
- if (!fEMCalCells)
+ if (!fVCells)
return;
- Int_t ncells = fEMCalCells->GetNumberOfCells();
- Int_t mcel = 0;
+ Int_t ncells = fVCells->GetNumberOfCells();
+ Int_t mcel = 0;//, maxcelid=-1;
+ Double_t maxcellE = 0;//, maxcellEta=0, maxcellPhi=0;
for(Int_t icell = 0; icell<ncells; icell++){
- Int_t absID = TMath::Abs(fEMCalCells->GetCellNumber(icell));
+ Int_t absID = TMath::Abs(fVCells->GetCellNumber(icell));
AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
Float_t eta=-1, phi=-1;
if(!fGeom){
if(!fGeom)
return;
/*if(!fIsMC)*/fGeom->EtaPhiFromIndex(absID,eta,phi);
+ if(maxcellE<fVCells->GetCellAmplitude(absID)){
+ maxcellE = fVCells->GetCellAmplitude(absID);
+ /*maxcellEta = eta;
+ maxcellPhi = phi;
+ maxcelid = absID;*/
+ }
Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
mycell->fAbsID = absID;
- mycell->fE = fEMCalCells->GetCellAmplitude(absID);
- mycell->fEt = fEMCalCells->GetCellAmplitude(absID)*TMath::Sin(theta);
+ mycell->fE = fVCells->GetCellAmplitude(absID);
+ mycell->fEt = fVCells->GetCellAmplitude(absID)*TMath::Sin(theta);
mycell->fEta = eta;
mycell->fPhi = phi;
- mycell->fTime = fEMCalCells->GetCellTime(absID);
+ mycell->fTime = fVCells->GetCellTime(absID);
}
+ if(this->fDebug)
+ printf("::FillMyCells() returning...\n\n");
}
//________________________________________________________________________
void AliAnalysisTaskEMCALPhoton::FillMyClusters()
{
// Fill clusters.
+ if(this->fDebug)
+ printf("::FillMyClusters() starting\n");
- if (!fCaloClusters)
+ if (!fCaloClusters){
+ printf("CaloClusters is empty!\n");
return;
+ }
Int_t nclus = fCaloClusters->GetEntries();
- Int_t mcl = 0;
+ if(0==nclus)
+ printf("CaloClusters has ZERO entries\n");
+ Int_t mcl = 0, maxcelid=-1;
+ Double_t maxcellE=0, maxcellEtac=0,maxcellPhic=0;
for(Int_t ic=0; ic < nclus; ic++){
- AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
+ AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(ic));
if(!clus)
continue;
if(!clus->IsEMCAL())
continue;
if(clus->E() < fClusThresh)
continue;
+ if(fDebug)
+ printf("cluster %d survived\n", ic);
Float_t pos[3];
clus->GetPosition(pos);
TVector3 cpos(pos);
Short_t id = -1;
myclus->fEmax = GetMaxCellEnergy( clus, id);
myclus->fIdmax = id;
- myclus->fTmax = fEMCalCells->GetCellTime(id);
+ if(maxcellE < myclus->fEmax){
+ maxcellE = myclus->fEmax;
+ maxcelid = id;
+ maxcellEtac = cpos.Eta();
+ maxcellPhic = cpos.Phi();
+ }
+ myclus->fTmax = fVCells->GetCellTime(id);
myclus->fEcross = GetCrossEnergy( clus, id);
myclus->fDisp = clus->GetDispersion();
myclus->fM20 = clus->GetM20();
myclus->fCellsAbsId = cellsAbsID;
myclus->fMcLabel = clus->GetLabel();
Int_t matchIndex = clus->GetTrackMatchedIndex();
- if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
- myclus->fTrEp = -1;
- continue;
- }
- AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
- if(!track){
+ if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
myclus->fTrEp = -1;
continue;
}
- if(!fPrTrCuts){
- myclus->fTrEp = -1;
- continue;
- }
- if(!fPrTrCuts->IsSelected(track)){
- myclus->fTrEp = -1;
+ AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
+ if(!track)
continue;
+ if(fESD){
+ if(!fPrTrCuts)
+ continue;
+ if(!fPrTrCuts->IsSelected(track))
+ continue;
}
+
myclus->fTrEp = clus->E()/track->P();
+ myclus->fTrDedx = track->GetTPCsignal();
+ }
+ if(this->fDebug){
+ printf(" ---===+++ Max Cell among clusters: id=%d, E=%1.2f, eta-clus=%1.2f, phi-clus=%1.2f\n",maxcelid,maxcellE,maxcellEtac,maxcellPhic);
+ printf("::FillMyClusters() returning...\n\n");
}
}
void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
{
// Fill clusters.
+ if(this->fDebug)
+ printf("::FillMyAltClusters() starting\n");
if(!fCaloClustersNew)
return;
Int_t nclus = fCaloClustersNew->GetEntries();
Int_t mcl = 0;
for(Int_t ic=0; ic < nclus; ic++){
- AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClustersNew->At(ic));
+ AliVCluster *clus = static_cast<AliVCluster*>(fCaloClustersNew->At(ic));
if(!clus)
continue;
if(!clus->IsEMCAL())
myclus->fCellsAbsId = cellsAbsID;
myclus->fMcLabel = clus->GetLabel();
Int_t matchIndex = clus->GetTrackMatchedIndex();
- if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
+ if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
myclus->fTrEp = -1;
continue;
}
- AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
+ AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
if(!track){
myclus->fTrEp = -1;
continue;
}
myclus->fTrEp = clus->E()/track->P();
}
+ if(this->fDebug)
+ printf("::FillMyAltClusters() returning...\n\n");
}
//________________________________________________________________________
void AliAnalysisTaskEMCALPhoton::FillIsoTracks()
{
// Fill high pt tracks.
+ if(this->fDebug)
+ printf("::FillIsoTracks() starting\n");
if(!fSelPrimTracks)
return;
Int_t ntracks = fSelPrimTracks->GetEntries();
Int_t imtr = 0;
for(Int_t it=0;it<ntracks; it++){
- AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(it));
+ AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(it));
if(!track)
continue;
- /*if(track->Pt()<3)
+ /*if(track->Phi()<1.0 || track->Phi()>3.55)
continue;*/
- AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
- if(track->Phi()<1.0 || track->Phi()>3.55)
- continue;
if(TMath::Abs(track->Eta())>1.1)
continue;
+ /*if(track->Pt()<3)
+ continue;*/
+ AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
mtr->fPt = track->Pt();
mtr->fEta = track->Eta();
mtr->fPhi = track->Phi();
mtr->fDedx = track->GetTPCsignal();
mtr->fMcLabel = track->GetLabel();
}
+ if(this->fDebug)
+ printf("::FillIsoTracks() returning...\n\n");
}
//________________________________________________________________________
void AliAnalysisTaskEMCALPhoton::GetMcParts()
{
- // Get MC particles.
+ // Get MC particles.
+ if(this->fDebug)
+ printf("::GetMcParts() starting\n");
- if (!fStack)
+ if (!this->fStack && !fAODMCParticles)
return;
-
- const AliVVertex *evtVtx = fMCEvent->GetPrimaryVertex();
- if (!evtVtx)
- return;
- Int_t ipart = 0;
- Int_t nTracks = fStack->GetNtrack();
+ if(this->fDebug)
+ printf("either stack or aodmcpaticles exists\n");
+ const AliVVertex *evtVtx = 0;
+ if(this->fStack)
+ evtVtx = fMCEvent->GetPrimaryVertex();
+ else
+ printf("no such thing as mc vvtx\n");
+ /*if (!evtVtx)
+ return;*/
+ if(this->fDebug)
+ printf("mc vvertex ok\n");
+ Int_t nTracks = 0;
+ if(this->fStack)
+ nTracks = this->fStack->GetNtrack();
+ else if(fAODMCParticles)
+ nTracks = fAODMCParticles->GetEntriesFast();
+ TParticle *mcP = 0;
+ AliAODMCParticle *amcP = 0;
+ if(this->fDebug)
+ printf("loop in the %d mc particles starting\n",nTracks);
for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
- TParticle *mcP = static_cast<TParticle*>(fStack->Particle(iTrack));
- if (!mcP)
- continue;
+ if(this->fStack)
+ mcP = dynamic_cast<TParticle*>(this->fStack->Particle(iTrack));
+ if(fAODMCParticles)
+ amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
+
// primary particle
- Double_t dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
- (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
- (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
- if(dR > 0.5)
+ Double_t dR = 0;
+ if(mcP){
+ dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
+ (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
+ (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
+ }
+ if((dR > 0.5))
continue;
+
// kinematic cuts
- Double_t pt = mcP->Pt() ;
+ Double_t pt = 0;
+ Double_t eta = 0;
+ Double_t phi = 0;
+ Int_t mother = -1;
+ Int_t pdgcode = 0;
+ if(mcP){
+ pt = mcP->Pt() ;
+ eta = mcP->Eta();
+ phi = mcP->Phi();
+ mother = mcP->GetMother(0);
+ pdgcode = mcP->GetPdgCode();
+ } else {
+ if(amcP){
+ pt = amcP->Pt();
+ eta = amcP->Eta();
+ phi = amcP->Phi();
+ mother = amcP->GetMother();
+ pdgcode = amcP->GetPdgCode();
+ } else
+ continue;
+ }
if (pt<0.5)
continue;
- Double_t eta = mcP->Eta();
- if (eta<-0.7||eta>0.7)
+
+ if (TMath::Abs(eta)>0.7)
continue;
- Double_t phi = mcP->Phi();
+
if (phi<1.0||phi>3.3)
continue;
+
+ if (mother!=6 && mother!=7 )
+ continue;
+
+
// pion or eta meson or direct photon
- if(mcP->GetPdgCode() == 111) {
- } else if(mcP->GetPdgCode() == 221) {
- } else if(mcP->GetPdgCode() == 22 ) {
- } else
- continue;
- bool checkIfAlreadySaved = false;
- for(Int_t imy=0;imy<fMyMcParts->GetEntries();imy++){
- AliPhotonMcPartObj *mymc = static_cast<AliPhotonMcPartObj*>(fMyMcParts->At(imy));
- if(!mymc)
- continue;
- if(mymc->fLabel == iTrack)
- checkIfAlreadySaved = true;
- }
- if(!checkIfAlreadySaved)
- FillMcPart(mcP, ipart++, iTrack);
- for(Int_t id=mcP->GetFirstDaughter(); id <= mcP->GetLastDaughter(); id++){
- if(id<=mcP->GetMother(0))
- continue;
- if(id<0 || id>nTracks)
- continue;
- TParticle *mcD = static_cast<TParticle*>(fStack->Particle(id));
- if(!mcD)
- continue;
- FillMcPart(mcD, ipart++, id);
+ if(pdgcode == 111) {
+ } else if(pdgcode == 221) {
+ } else if(pdgcode == 22 ) {
+ } else {
+ continue;
}
+
+ FillMcPart( fMyMcIndex++, iTrack);
}
+ if(this->fDebug)
+ printf("::GetMcParts() returning...\n\n");
}
//________________________________________________________________________
-void AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t ipart, Int_t itrack)
+void AliAnalysisTaskEMCALPhoton::FillMcPart( Int_t itrack, Int_t label)
{
// Fill MC particles.
-
- if(!mcP)
- return;
- TVector3 vmcv(mcP->Vx(),mcP->Vy(), mcP->Vz());
- AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(ipart));
- mcp->fLabel = itrack ;
- mcp->fPdg = mcP->GetPdgCode() ;
- mcp->fPt = mcP->Pt() ;
- mcp->fEta = mcP->Eta() ;
- mcp->fPhi = mcP->Phi() ;
+ Int_t nTracks = 0;
+ TParticle *mcP = 0;
+ AliAODMCParticle *amcP= 0;
+ if(this->fStack){
+ nTracks = this->fStack->GetNtrack();
+ mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
+ }
+ else if(fAODMCParticles){
+ nTracks = fAODMCParticles->GetEntriesFast();
+ amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(itrack));
+ }
+ if(this->fDebug)
+ printf("\t::FillMcParts() starting with label %d\n", itrack);
+ TVector3 vmcv;
+ if(mcP)
+ vmcv.SetXYZ(mcP->Vx(),mcP->Vy(), mcP->Vz());
+ else {
+ if(amcP)
+ vmcv.SetXYZ(amcP->Xv(),amcP->Yv(), amcP->Zv());
+ else
+ return;
+ }
+
+ AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(itrack));
+ mcp->fLabel = label ;
+ if(mcP){
+ mcp->fPdg = mcP->GetPdgCode() ;
+ mcp->fPt = mcP->Pt() ;
+ mcp->fEta = mcP->Eta() ;
+ mcp->fPhi = mcP->Phi() ;
+ mcp->fMother = mcP->GetMother(0) ;
+ mcp->fFirstD = mcP->GetFirstDaughter() ;
+ mcp->fLastD = mcP->GetLastDaughter() ;
+ mcp->fStatus = mcP->GetStatusCode();
+ }
+ if(amcP){
+ mcp->fPdg = amcP->GetPdgCode() ;
+ mcp->fPt = amcP->Pt() ;
+ mcp->fEta = amcP->Eta() ;
+ mcp->fPhi = amcP->Phi() ;
+ mcp->fMother = amcP->GetMother() ;
+ mcp->fFirstD = amcP->GetDaughter(0) ;
+ mcp->fLastD = amcP->GetDaughter(amcP->GetNDaughters()-1) ;
+ mcp->fStatus = amcP->GetStatus();
+ }
mcp->fVR = vmcv.Perp();
if(vmcv.Perp()>0){
mcp->fVEta = vmcv.Eta() ;
mcp->fVPhi = vmcv.Phi() ;
}
- mcp->fMother = mcP->GetMother(0) ;
+ if(itrack == 8){
+ mcp->fIso = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.4 , 0.2);
+ mcp->fIso3 = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.3 , 0.2);
+ }
+ if(this->fDebug)
+ printf("\t::FillMcParts(): label=%d, pdg=%d, pt=%1.1f, mom=%d, 1stD=%d,last=%d\n\t::FillMcParts() returning...\n\n", mcp->fLabel,mcp->fPdg,mcp->fPt,mcp->fMother,mcp->fFirstD,mcp->fLastD);
+ for(Int_t id=mcp->fFirstD; id < mcp->fLastD; id++){
+ if(id<=mcp->fMother)
+ continue;
+ if(id<0 || id>nTracks)
+ continue;
+ FillMcPart( fMyMcIndex++, id);
+ }
+
}
+//________________________________________________________________________
+Double_t AliAnalysisTaskEMCALPhoton::GetMcIsolation( Int_t itrack, Double_t radius, Double_t ptcut) const
+{
+ if(this->fDebug){
+ printf("\t\t::GetMcIsolation() starting\n");
+ //printf("\t\t incoming particle: PDG = %d, itrack=%d;\n",mcP->GetPdgCode(),itrack);
+ }
+ if (!this->fStack && !this->fAODMCParticles && this->fDebug){
+ printf("\t\t\tNo MC stack/array!\n");
+ return -1;
+ }
+ if(itrack<6 || itrack>8)
+ return -1;
+ if(this->fDebug)
+ printf("\t\t\tparticle of interest selected\n");
+ TParticle *mcP = 0;
+ AliAODMCParticle *amcP = 0;
+ Int_t pdgcode = 0;
+ Float_t eta = 0;
+ Float_t phi = 0;
+ Double_t sumpt=0;
+ Float_t dR;
+ Int_t nparts = 0;
+ if(this->fStack){
+ nparts = fStack->GetNtrack();
+ mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
+ eta = mcP->Eta();
+ phi = mcP->Phi();
+ pdgcode = mcP->GetPdgCode();
+ }
+ if(this->fAODMCParticles){
+ nparts = fAODMCParticles->GetEntriesFast();
+ amcP = dynamic_cast<AliAODMCParticle*>(this->fAODMCParticles->At(itrack));
+ if(amcP){
+ eta = amcP->Eta();
+ phi = amcP->Phi();
+ pdgcode = amcP->GetPdgCode();
+ }
+ }
+ if(pdgcode!=22)
+ return -1;
+ TParticle *mcisop = 0;
+ AliAODMCParticle *amcisop = 0;
+ for(Int_t ip = 0; ip<nparts; ip++){
+ if(ip==itrack)
+ continue;
+ if(this->fStack)
+ mcisop = dynamic_cast<TParticle*>(this->fStack->Particle(ip));
+ if(fAODMCParticles)
+ amcisop = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(ip));
+ Int_t status = 0;
+ Int_t mother = 0;
+ Float_t pt = 0;
+ Float_t isophi = -99;
+ Float_t isoeta = -99;
+ TVector3 vmcv;
+ if(mcisop){
+ status = mcisop->GetStatusCode();
+ mother = mcisop->GetMother(0);
+ pt = mcisop->Pt();
+ isophi = mcisop->Phi();
+ isoeta = mcisop->Eta();
+ vmcv.SetXYZ(mcisop->Vx(),mcisop->Vy(), mcisop->Vz());
+ }
+ else {
+ if(amcisop){
+ status = amcisop->GetStatus();
+ mother = amcisop->GetMother();
+ pt = amcisop->Pt();
+ isophi = amcisop->Phi();
+ isoeta = amcisop->Eta();
+ vmcv.SetXYZ(amcisop->Xv(),amcisop->Yv(), amcisop->Zv());
+ }
+ else
+ continue;
+ }
+ if(status!=1)
+ continue;
+ if(mother == itrack)
+ continue;
+ if(pt<ptcut)
+ continue;
+ if(vmcv.Perp()>1)
+ continue;
+ dR = TMath::Sqrt((phi-isophi)*(phi-isophi)+(eta-isoeta)*(eta-isoeta));
+ if(dR>radius)
+ continue;
+ sumpt += pt;
+ }
+ if(this->fDebug)
+ printf("\t\t::GetMcIsolation() returning value %f ...\n\n",sumpt);
+ return sumpt;
+ }
//________________________________________________________________________
Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
{
// Compute isolation based on tracks.
-
+ if(this->fDebug)
+ printf("\t::GetTrackIsolation() starting\n");
+
Double_t trkIsolation = 0;
Double_t rad2 = radius*radius;
Int_t ntrks = fSelPrimTracks->GetEntries();
continue;
trkIsolation += track->Pt();
}
+ if(this->fDebug)
+ printf("\t::GetTrackIsolation() returning\n\n");
return trkIsolation;
}
Double_t lowphi = phi - R;
Double_t et = 0;
for(Int_t itr=0; itr<ntracks; itr++){
- AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(itr));
+ AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(itr));
if(!track)
continue;
if(track->Pt()<minpt)
Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
{
// Calculate the energy of cross cells around the leading cell.
-
- AliVCaloCells *cells = 0;
- cells = fESD->GetEMCALCells();
- if (!cells)
+ if(!fVCells)
return 0;
if (!fGeom)
continue;
if ( (aphidiff==1 && aetadiff==0) ||
(aphidiff==0 && aetadiff==1) ) {
- crossEnergy += cells->GetCellAmplitude(cellAbsId);
+ crossEnergy += fVCells->GetCellAmplitude(cellAbsId);
}
}
// Get maximum energy of attached cell.
id = -1;
-
- AliVCaloCells *cells = 0;
- cells = fESD->GetEMCALCells();
- if (!cells)
+ if(!fVCells)
return 0;
Double_t maxe = 0;
Int_t ncells = cluster->GetNCells();
for (Int_t i=0; i<ncells; i++) {
- Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
+ Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
if (e>maxe) {
maxe = e;
id = cluster->GetCellAbsId(i);
/* if(fIsGrid)
return;*/
if (fTree) {
- TFile *f = OpenFile(1);
+ printf("***tree %s being saved***\n",fTree->GetName());
+ TFile *f = OpenFile(2);
TDirectory::TContext context(f);
if (f)
fTree->Write();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff