+//________________________________________________________________
+void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
+ //Estimate background rho by means of integrating track pT outside identified jet cones
+
+ rhoMedian = 0;
+
+ //identify leading jet in the full track acceptance
+ Int_t idxLeadingJet = -1;
+ Double_t pTleading = 0.0;
+ AliAODJet* jet = NULL;
+
+ for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
+ jet = (AliAODJet*)(listJet->At(ij));
+ if(!jet) continue;
+ if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
+ if(jet->Pt() > pTleading){
+ idxLeadingJet = ij;
+ pTleading = jet->Pt();
+ }
+ }
+
+ Int_t njets = 0;
+ static Double_t jphi[100];
+ static Double_t jeta[100];
+ static Double_t jRsquared[100];
+
+ for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
+
+ jet = (AliAODJet*)(listJet->At(ij));
+ if(!jet) continue;
+ if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
+
+ if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
+ jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
+ jeta[njets] = jet->Eta();
+ jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
+ njets++;
+ }
+ }
+
+
+ static Double_t nOutCone[10][4];
+ static Double_t sumPtOutOfCone[10][4];
+
+ for(Int_t ie=0; ie<fnEta; ie++){
+ for(Int_t ip=0; ip<fnPhi; ip++){
+ nOutCone[ip][ie] = 0.0; //initialize counter
+ sumPtOutOfCone[ip][ie] = 0.0;
+ }
+ }
+
+ Double_t rndphi, rndeta;
+ Double_t rndphishift, rndetashift;
+ Double_t dphi, deta;
+ Bool_t bIsInCone;
+
+ for(Int_t it=0; it<fnTrials; it++){
+
+ rndphi = fRandom->Uniform(0, fPhiSize);
+ rndeta = fRandom->Uniform(0, fEtaSize);
+
+ for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
+ rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
+ for(Int_t ie=0; ie<fnEta; ie++){
+ rndetashift = rndeta + ie*fEtaSize - fEtaSize;
+
+ bIsInCone = 0; //tag if trial is in the jet cone
+ for(Int_t ij=0; ij<njets; ij++){
+ deta = jeta[ij] - rndetashift;
+ dphi = RelativePhi(rndphishift,jphi[ij]);
+ if((dphi*dphi + deta*deta) < jRsquared[ij]){
+ bIsInCone = 1;
+ break;
+ }
+ }
+ if(!bIsInCone) nOutCone[ip][ie]++;
+ }
+ }
+ }
+
+
+ //Sum particle pT outside jets in cells
+ Int_t npart = listPart->GetEntries();
+ Int_t phicell,etacell;
+ AliVParticle *part = NULL;
+ for(Int_t ip=0; ip < npart; ip++){
+
+ part = (AliVParticle*) listPart->At(ip);
+ if(!part) continue;
+ if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
+
+ bIsInCone = 0; //init
+ for(Int_t ij=0; ij<njets; ij++){
+ dphi = RelativePhi(jphi[ij], part->Phi());
+ deta = jeta[ij] - part->Eta();
+ if((dphi*dphi + deta*deta) < jRsquared[ij]){
+ bIsInCone = 1;
+ break;
+ }
+ }
+ if(!bIsInCone){
+ phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
+ etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
+ sumPtOutOfCone[phicell][etacell]+= part->Pt();
+ }
+ }
+
+ static Double_t rhoInCells[20];
+ Double_t relativeArea;
+ Int_t nCells=0;
+ Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
+ for(Int_t ip=0; ip<fnPhi; ip++){
+ for(Int_t ie=0; ie<fnEta; ie++){
+ relativeArea = nOutCone[ip][ie]/fnTrials;
+ //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
+
+ bufferArea += relativeArea;
+ bufferPt += sumPtOutOfCone[ip][ie];
+ if(bufferArea > fJetFreeAreaFrac){
+ rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
+
+ if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
+ else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
+
+ bufferArea = 0.0;
+ bufferPt = 0.0;
+ nCells++;
+ }
+ }
+ }
+
+
+ if(nCells>0){
+ rhoMedian = TMath::Median(nCells, rhoInCells);
+ if(mode==1){ //mc data
+ fhEntriesToMedianGen->Fill(nCells);
+ }else{ //real data
+ fhEntriesToMedian->Fill(nCells);
+ }
+ }
+
+}
+
+//__________________________________________________________________
+void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
+
+ rhoPerpCone = 0.0; //init
+
+ if(!listJet) return;
+ Int_t njet = listJet->GetEntries();
+ Int_t npart = listPart->GetEntries();
+ Double_t pTleading = 0.0;
+ Double_t phiLeading = 1000.;
+ Double_t etaLeading = 1000.;
+
+ for(Int_t jsig=0; jsig < njet; jsig++){
+ AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
+ if(!signaljet) continue;
+ if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
+ if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
+ pTleading = signaljet->Pt();
+ phiLeading = signaljet->Phi();
+ etaLeading = signaljet->Eta();
+ }
+ }
+
+ Double_t phileftcone = phiLeading + TMath::Pi()/2;
+ Double_t phirightcone = phiLeading - TMath::Pi()/2;
+ Double_t dp, de;
+
+ for(Int_t ip=0; ip < npart; ip++){
+
+ AliVParticle *part = (AliVParticle*) listPart->At(ip);
+ if(!part){
+ continue;
+ }
+
+ dp = RelativePhi(phileftcone, part->Phi());
+ de = etaLeading - part->Eta();
+ if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
+
+ dp = RelativePhi(phirightcone, part->Phi());
+ if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
+
+ }