+//________________________________________________________________________
+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;
+ }