+
+//___________________________________________________________________________________________________________________
+void AliEmcalJetTask::FillJetConstituents(std::vector<fastjet::PseudoJet>& constituents,AliEmcalJet *jet,Double_t vertex[3],Int_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc){
+ nt = 0;
+ nc = 0;
+ neutralE = 0;
+ maxCh = 0;
+ maxNe = 0;
+ gall = 0;
+ gemc =0;
+ cemc = 0;
+ ncharged = 0;
+ nneutral = 0;
+ mcpt = 0;
+ emcpt = 0;
+ AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
+ for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
+ Int_t uid = constituents[ic].user_index();
+ AliDebug(3,Form("Processing constituent %d", uid));
+ if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
+ ++gall;
+ Double_t gphi = constituents[ic].phi();
+ if (gphi<0)
+ gphi += TMath::TwoPi();
+ gphi *= TMath::RadToDeg();
+ Double_t geta = constituents[ic].eta();
+ if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
+ (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
+ ++gemc;
+ } else if ((uid > 0) && fTracks) { // track constituent
+ Int_t tid = uid - 100;
+ AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
+ if (!t) {
+ AliError(Form("Could not find track %d",tid));
+ continue;
+ }
+ if (jetCount < fMarkConst) {
+ AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
+ t->SetBit(fJetType);
+ }
+ Double_t cEta = t->Eta();
+ Double_t cPhi = t->Phi();
+ Double_t cPt = t->Pt();
+ Double_t cP = t->P();
+ if (t->Charge() == 0) {
+ neutralE += cP;
+ ++nneutral;
+ if (cPt > maxNe)
+ maxNe = cPt;
+ } else {
+ ++ncharged;
+ if (cPt > maxCh)
+ maxCh = cPt;
+ }
+ if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
+ mcpt += cPt;
+
+ if (cPhi<0)
+ cPhi += TMath::TwoPi();
+ cPhi *= TMath::RadToDeg();
+ if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
+ (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
+ emcpt += cPt;
+ ++cemc;
+ }
+
+ jet->AddTrackAt(tid, nt);
+ ++nt;
+ } else if (fClus) { // cluster constituent
+ Int_t cid = -uid - 100;
+ AliVCluster *c = 0;
+ Double_t cEta=0,cPhi=0,cPt=0,cP=0;
+ if (fIsEmcPart) {
+ AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
+ if (!ep)
+ continue;
+ c = ep->GetCluster();
+ if (!c)
+ continue;
+ if (jetCount < fMarkConst)
+ ep->SetBit(fJetType);
+ cEta = ep->Eta();
+ cPhi = ep->Phi();
+ cPt = ep->Pt();
+ cP = ep->P();
+ } else {
+ c = static_cast<AliVCluster*>(fClus->At(cid));
+ if (!c)
+ continue;
+ if (jetCount < fMarkConst)
+ c->SetBit(fJetType);
+ TLorentzVector nP;
+ c->GetMomentum(nP, vertex);
+ cEta = nP.Eta();
+ cPhi = nP.Phi();
+ cPt = nP.Pt();
+ cP = nP.P();
+ }
+
+ neutralE += cP;
+ if (cPt > maxNe)
+ maxNe = cPt;
+
+ if (c->GetLabel() > fMinMCLabel) // MC particle
+ mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
+
+ if (cPhi<0)
+ cPhi += TMath::TwoPi();
+ cPhi *= TMath::RadToDeg();
+ if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
+ (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
+ emcpt += cPt;
+ ++cemc;
+ }
+
+ jet->AddClusterAt(cid, nc);
+ ++nc;
+ ++nneutral;
+ } else {
+ AliError(Form("%s: No logical way to end up here.", GetName()));
+ continue;
+ }
+
+
+
+
+ }
+
+
+}
+//______________________________________________________________________________________