Index: trunk/code/CMS_2014_I1299142.cc
===================================================================
--- trunk/code/CMS_2014_I1299142.cc	(revision 475)
+++ trunk/code/CMS_2014_I1299142.cc	(revision 476)
@@ -1,602 +1,602 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/ChargedFinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
 
 namespace Rivet {
 
 using namespace fastjet;
 	
 	struct FourMomComp {
         bool operator() (const FourMomentum& p1, const FourMomentum& p2) const
        {return p1.pT()<p2.pT();}
     };
 
 
     struct jetcomp{
       bool operator() (const Jet jet1, const Jet jet2) const
       {
         if (jet1.momentum().pT() < jet2.momentum().pT()) return true;
         else return false;
       }
     };
 
     struct My4Mom {
 	    double pt;
 	    double mdelta;
 	    double phi;
 	    double y;
     } ;
 
     struct MyPart {
 	    double pt;
 	    double mdelta;
 	    double phi;
 	    double y;
 	    PdgId  id;
     } ;
 
     struct MyDist {
 	    double dR;
 	    size_t ipart;
 	    size_t ighost;
     } ;
     
     bool MyDistComp (const MyDist& dist1, const MyDist& dist2) {
 	    return dist1.dR < dist2.dR; 
     }
 
 
 /*#if MODE==0
   class CMS_2014_I1299142_4MOMSUB : public Analysis {
 #elif MODE==1
   class CMS_2014_I1299142_GRIDSUB2 : public Analysis {
 #elif MODE==3
   class CMS_2014_I1299142_CONSTSUB : public Analysis {
 #endif*/
   class CMS_2014_I1299142 : public Analysis {
   public:
 
     /// Constructor
 /*#if MODE==0
     CMS_2014_I1299142_4MOMSUB() : Analysis("CMS_2014_I1299142_4MOMSUB")
 #elif MODE==1
     CMS_2014_I1299142_GRIDSUB2() : Analysis("CMS_2014_I1299142_GRIDSUB2")
 #elif MODE==3
     CMS_2014_I1299142_CONSTSUB() : Analysis("CMS_2014_I1299142_CONSTSUB")
 #endif*/
     CMS_2014_I1299142() : Analysis("CMS_2014_I1299142")
     {
       _jetR=0.3; _trackptmin=1.; _etamin=0.3; _etamax=2.; _ptmin=100.; _ptmax=300.; _Nptbins=3;
-	_etaloose=2.5; _ptminloose=70.; _dRloose=0.6;
+	_etaloose=2.3; _ptminloose=90.; _dRloose=1.0;
       _nphibins = 120; _netabins = 160;
     }
 
 
     /// Book histograms and initialise projections before the run
     void init() {
 
-	 _ptedges += 100.0, 120.0, 150.0, 300.0;
+	 _ptedges = {100.0, 120.0, 150.0, 300.0};
       _Njets.resize(_Nptbins+1, 0.0);
 
 	 FinalState fs(-_etaloose, _etaloose, 0.*GeV);
        addProjection(fs, "FS");
        ChargedFinalState cfs(-5.0, 5.0, _trackptmin*GeV);
        addProjection(cfs, "CFS");
 	
        FastJets fj(fs, FastJets::ANTIKT, _jetR);
        addProjection(fj, "Jets");
 
 	 for (size_t i = 0; i < _Nptbins+1; ++i) {
 		 _h_xi_pp[i]    = bookHisto1D(i+1, 1, 1);
 		 _h_xi_PbPb[i]  = bookHisto1D(i+1, 2, 1);
 		 _h_xi_ratio[i] = bookHisto1D(i+1, 3, 1);
 		 _h_pt_pp[i]    = bookHisto1D(i+1, 1, 2);
 		 _h_pt_PbPb[i]  = bookHisto1D(i+1, 2, 2);
 		 _h_pt_diff[i]  = bookHisto1D(i+1, 3, 2);
 	 }
 	
     }
 
 
 
     double ptsmear(const double pttrue, const double cent) {
       double c(0.0246);
       double s(1.213);
       double n;
       if (cent < 0.  || cent > 0.5) n=0.001;
       if (cent > 0.3 && cent < 0.5) n=3.88;
       if (cent > 0.1 && cent < 0.3) n=5.10;
       if (cent > 0.  && cent < 0.1) n=5.23;
       n=5.23;
       double sigma(sqrt(c*c+s*s/pttrue+n*n/(pttrue*pttrue)));
       double r1(1.0*rand()/RAND_MAX), r2(1.0*rand()/RAND_MAX);
       double fac(max(sqrt(-2.*log(r1))*cos(2.*M_PI*r2)*sigma+1.,0.));
       return fac*pttrue;
     }    
 
     
     /// extract thermal momenta, that should be subtracted, from HepMC event
     vector<FourMomentum> extractThermalMomenta(const Event & event) {
 	 vector<FourMomentum> thermom;
         foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
             FourMomentum mom(p->momentum());
                if (p->status() == 3) {
                    thermom.push_back(mom);
                }
 	 }
 	 return thermom;
     }
     
     
     /// build map from dummy particles to thermal momenta to speed up subtraction
     map<FourMomentum, FourMomentum, FourMomComp> buildThMomMap(const Event & event, const vector<FourMomentum> * thermom) {
 	    map<FourMomentum, FourMomentum, FourMomComp> thmap;
 	    size_t j(0);
            foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
                   FourMomentum mom(p->momentum());
 			if (fuzzyEquals(mom.E(), 1e-6)) {
 				if (j == thermom->size()) {
 					cout<<"Error: number of dummies does not match number of subtractions, will veto Event"<<endl;
 					thmap.clear();
 					return thmap;
 				}
 				FourMomentum sc((*thermom)[j]);
 				thmap[mom]=sc;
 				j++;
 			}
 	    }
 	    return thmap;
     }
    
     
     /// 4-momentum subtraction
     FourMomentum SubtractJetMom(Jet jet, map<FourMomentum, FourMomentum, FourMomComp> * tmmap) {
         if (tmmap->empty()) return jet.momentum();
         FourMomentum sub(0.,0.,0.,0.);
         foreach (Particle part, jet.constituents()) {
             if (fuzzyEquals(part.E(), 1e-6)) {	 
                     map<FourMomentum, FourMomentum, FourMomComp>::iterator mapit;
                     mapit = (*tmmap).find(part.momentum());
                     if (mapit == (*tmmap).end()) cout<<"Error: did not find matching scattering centre in map.\n"<<part.momentum()<<endl;
                     else sub += mapit->second;
 	     }
         }
         FourMomentum jetmom(jet.momentum());
         //cout<<"Original momentum: "<<jetmom<<endl;
         jetmom -= sub;
         //cout<<"subtracted momentum: "<<jetmom<<endl;
         //cout<<"Subtracted pt: "<<subpt<<"    "<<pjet.pt()<<endl;
         return jetmom;
     }
 
     
 /// initialise and fill grid for grid subtration    
     vector<vector<FourMomentum> >  fillGrid(const Event & event, std::set<std::pair<size_t,size_t> > & gridevent) {	    
 	    // initialise grid
            vector<vector<FourMomentum> >  grid;
            grid.resize(_nphibins);
            for (size_t i = 0; i < _nphibins; ++i) {
                   grid[i].resize(_netabins);
 		    for (size_t k = 0; k < _nphibins; ++k) {
 			    grid[i][k] = FourMomentum(0., 0., 0., 0.);
 		    }
            }
 
       // fill grid
         foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
             FourMomentum mom(p->momentum());
             if (fabs(mom.eta()) < 4.) {
                size_t phibin, etabin;
                phibin = int(mom.phi(ZERO_2PI)*_nphibins/(2.*M_PI));
                etabin = int((mom.eta() + 4.)*_netabins/8.);
                if (phibin < 0 || phibin > _nphibins-1) {std::cout<<"Error: "<<mom.phi(ZERO_2PI)<<"  --> "<<phibin<<std::endl; exit(1);}
                if (etabin < 0 || etabin > _netabins-1) {std::cout<<"Error: "<<mom.eta()<<"  --> "<<etabin<<std::endl; exit(1);}
                if (p->status() == 3) {
                    grid[phibin][etabin].setE(grid[phibin][etabin].E() - mom.E());
                    grid[phibin][etabin].setPx(grid[phibin][etabin].px() - mom.px());
                    grid[phibin][etabin].setPy(grid[phibin][etabin].py() - mom.py());
                    grid[phibin][etabin].setPz(grid[phibin][etabin].pz() - mom.pz());
                }
                else {
                    grid[phibin][etabin].setE(grid[phibin][etabin].E() + mom.E());
                    grid[phibin][etabin].setPx(grid[phibin][etabin].px() + mom.px());
                    grid[phibin][etabin].setPy(grid[phibin][etabin].py() + mom.py());
                    grid[phibin][etabin].setPz(grid[phibin][etabin].pz() + mom.pz());    
                }
                /*double cellphi((phibin+0.5)*2*M_PI/_nphibins), celleta((etabin+0.5)*8./_netabins-4.);
                double celltheta(2.*atan(exp(-celleta)));
                if (p->status() == 3) {
                    double rho(grid[phibin][etabin].E() - mom.E());
                    grid[phibin][etabin].setE(rho);
                    grid[phibin][etabin].setPx(fabs(rho)*sin(celltheta)*cos(cellphi));
                    grid[phibin][etabin].setPy(fabs(rho)*sin(celltheta)*sin(cellphi));
                    grid[phibin][etabin].setPz(fabs(rho)*cos(celltheta));
                }
                else {
 		     // use this for standard recoil treatment
                    double rho(grid[phibin][etabin].E() + mom.E());
                    grid[phibin][etabin].setE(rho);
                    grid[phibin][etabin].setPx(rho*sin(celltheta)*cos(cellphi));
                    grid[phibin][etabin].setPy(rho*sin(celltheta)*sin(cellphi));
                    grid[phibin][etabin].setPz(rho*cos(celltheta));
 		 }*/
                gridevent.insert(std::pair<size_t,size_t>(phibin,etabin));
 	     }
 	 }
 		return grid;
     }  
 
     /// event-wise Constituent subtraction based on Particles
     Particles ConstSubPartEvent(const Particles particles, PseudoJets pjghosts) {
 	    // sort constituents into two vectors: particles and ghosts (subtraction momenta)
 	    Particles subevent;
 	    vector<MyPart> parts;
 	    vector<My4Mom> ghosts;
 	    
 	    foreach (Particle p, particles) {
 			FourMomentum partmom = p.momentum();
 			if (fuzzyEquals(partmom.E(),1e-6)) continue;
 			if (partmom.E() - fabs(partmom.pz()) > 0.) {
 				MyPart part;
 				part.id = p.pid();
 				part.pt = partmom.pT();
 				part.mdelta = sqrt(partmom.mass2() + sqr(partmom.pT())) - partmom.pT();
 				part.phi = partmom.p3().phi();
 				part.y = partmom.rapidity();
 				//if (abs(part.y) > 5.) {
 				//	cout<<"particle rapidity: "<<part.y<<"     "<<partmom.E()-partmom.pz()<<endl;
 				//	cout<<partmom<<"   "<<part.pt<<endl;
 				//}
 				parts.push_back(part);
 			}
 			else {
 				subevent.push_back(p);
 			}
 	    }
 	    
 	    foreach (PseudoJet pj, pjghosts) {
 			FourMomentum submom = Get4Mom(pj);
 			if (submom.E() - fabs(submom.pz()) > 0.) {
 				My4Mom ghost;
 				ghost.pt = submom.pT();
 				ghost.mdelta = sqrt(submom.mass2() + sqr(submom.pT())) - submom.pT();
 				ghost.phi = submom.p3().phi();
 				ghost.y = submom.rapidity();
 				//cout<<"ghost rapidity: "<<ghost.y<<endl;
 				ghosts.push_back(ghost);
 			}
 			else {
 				PdgId ghostid;
 				if (rand()/RAND_MAX < 1./3.) ghostid = PID::PIPLUS;
 				else if (rand()/RAND_MAX < 2./3.) ghostid = PID::PIMINUS;
 				else ghostid = PID::PI0;
 				Particle pg(ghostid,Get4Mom(pj));
 				subevent.push_back(pg);
 			}
 	    }
 	    
         //cout<<"number of particles: "<<parts.size()<<endl;
         //cout<<"number of ghosts: "<<ghosts.size()<<endl;
 
 	  //create list with all particle-ghosts distances and sort it
 	  vector<MyDist> dists;
 	  for (size_t i=0; i < parts.size(); ++i) {
 		  for (size_t j=0; j < ghosts.size(); ++j) {
 			  //cout<<"i, j = "<<i<<" "<<j<<endl;
 			  MyDist dist;
 			  double Deltaphi(abs(parts[i].phi-ghosts[j].phi));
 			  if (Deltaphi > M_PI) Deltaphi = 2.*M_PI - Deltaphi;
 			  dist.dR = sqrt(sqr(Deltaphi) + sqr(parts[i].y-ghosts[j].y));
 			  dist.ipart = i;
 			  dist.ighost = j;
 			  //cout<<"distance: "<<dist.dR<<endl;
 			  dists.push_back(dist);
 		  }
 	  }
 	  //cout<<"number of pairs: "<<dists.size()<<endl;
 	  //dists.sort(MyDistComp);
 	  std::sort(dists.begin(),dists.end(),MyDistComp);
 	  
 	  //go through all particle-ghost pairs and re-distribute momentum and mass
 	  for (vector<MyDist>::iterator liter=dists.begin(); liter != dists.end(); ++liter) {
 		  //cout<<"dealing with dist "<<liter->dR<<endl;
 		  if (liter->dR > _dRloose) break;
 		  size_t pnum = liter->ipart;
 		  size_t gnum = liter->ighost;
 		  double ptp = parts[pnum].pt;
 		  double ptg = ghosts[gnum].pt;
 		  //cout<<"pts: "<<ptp<<" vs "<<ptg<<endl;
 		  if (ptp > ptg) {
 			  parts[pnum].pt -= ptg;
 			  ghosts[gnum].pt = 0.;
 		  }
 		  else {
 			  ghosts[gnum].pt -= ptp;
 			  parts[pnum].pt = 0.;
 		  }
 		  double mdp = parts[pnum].mdelta;
 		  double mdg = ghosts[gnum].mdelta;
 		  //cout<<"masses: "<<mdp<<" vs "<<mdg<<endl;
 		  if (mdp > mdg) {
 			  parts[pnum].mdelta -= mdg;
 			  ghosts[gnum].mdelta = 0.;
 		  }
 		  else {
 			  ghosts[gnum].mdelta -= mdp;
 			  parts[pnum].mdelta = 0.;
 		  }
 		  
 	  }
 	  
 	  //sum up resulting 4-momenta to get subtracted jet momentum
 	  int nparts(0);
 	  for (size_t i=0; i < parts.size(); ++i) {
 		  if (parts[i].pt > 0.) {
 			  FourMomentum mom((parts[i].pt+parts[i].mdelta)*cosh(parts[i].y), 
 						  parts[i].pt*cos(parts[i].phi),
 						  parts[i].pt*sin(parts[i].phi),
 						 (parts[i].pt+parts[i].mdelta)*sinh(parts[i].y));
 			  //cout<<"particle rapidity: "<<mom.rapidity()<<"  "<<parts[i].y<<endl;
 			  Particle outpart(parts[i].id,mom);
 			  subevent.push_back(outpart);
 			  nparts++;
 		  }
 	  }
 	  int nghost(0);
 	  for (size_t i=0; i < ghosts.size(); ++i) {
 		  if (ghosts[i].pt > 0.) {
 			  FourMomentum mom((ghosts[i].pt+ghosts[i].mdelta)*cosh(ghosts[i].y), 
 						  ghosts[i].pt*cos(ghosts[i].phi),
 						  ghosts[i].pt*sin(ghosts[i].phi),
 						 (ghosts[i].pt+ghosts[i].mdelta)*sinh(ghosts[i].y));
 			  //cout<<"ghost rapidity: "<<mom.rapidity()<<"   "<<ghosts[i].y<<endl;
 			  PdgId ghostid;
 			  if (rand()/RAND_MAX < 1./3.) ghostid = PID::PIPLUS;
 			  else if (rand()/RAND_MAX < 2./3.) ghostid = PID::PIMINUS;
 			  else ghostid = PID::PI0;
 			  Particle outpart(ghostid,mom);
 			  subevent.push_back(outpart);
 			  nghost++;	
 		}
 		  //else cout<<ghosts[i].mdelta<<endl;
 	  }
 	  //cout<<nghost<<" / "<<ghosts.size()+1<<" vs. "<<nparts<<" / "<<parts.size()<<endl;
         return subevent;
     }
 
 
 
 /// 4-momentum adapter
     FourMomentum Get4Mom(PseudoJet pjet){
         FourMomentum mom(pjet.E(), pjet.px(), pjet.py(), pjet.pz());
         return mom;
     }
     
     
     
 
     /// Perform the per-event analysis
     void analyze(const Event& event) {
       const double weight = event.weight();
       const double cent = (event.genEvent()->heavy_ion()?event.genEvent()->heavy_ion()->impact_parameter():-1.);
       if (cent > 0.1) vetoEvent;
       const ChargedFinalState cfs = applyProjection<ChargedFinalState>(event, "CFS");
       Particles tracks;
 	
       Jets jetssub;
 #if MODE==0      
 	tracks = cfs.particles();
       const Jets alljets = applyProjection<FastJets>(event, "Jets").jetsByPt(3.*GeV);
 
       vector<FourMomentum> thermom = extractThermalMomenta(event);
       map<FourMomentum, FourMomentum, FourMomComp> thmommap = buildThMomMap(event, &thermom);
       foreach (Jet jet, alljets) {
 	      FourMomentum jetsub = SubtractJetMom(jet, &thmommap);
 	      jetssub.push_back(Jet(jetsub));
       }
 #elif MODE==1
 	tracks = cfs.particles();
       PseudoJets pevent, pjets;
       std::set<std::pair<size_t,size_t> > gridevent;
       vector<vector<FourMomentum> >  grid = fillGrid(event, gridevent);
       for (set<pair<size_t,size_t> >::iterator it=gridevent.begin(); it!=gridevent.end(); ++it) {
           size_t i,k;
           i = it->first;
           k = it->second;
           if (grid[i][k].E() > 0. ) {
               PseudoJet part(grid[i][k].px(),grid[i][k].py(),grid[i][k].pz(),grid[i][k].E());
               pevent.push_back(part);
           }
       }
       JetDefinition jet_def(antikt_algorithm, _jetR);
       ClusterSequence cs(pevent, jet_def);
       pjets = sorted_by_pt(cs.inclusive_jets(3.*GeV));
       foreach (PseudoJet pjet, pjets) {
 	      jetssub.push_back(Jet(pjet));
       }
 #elif MODE==2
       const FinalState fs = applyProjection<FinalState>(event, "FS");
 	const Particles origevent = fs.particles();
       
 	PseudoJets ghosts;
 
 	foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
             FourMomentum mom(p->momentum());
 		if (fabs(mom.eta()) < _etaloose) {
 		      if (p->status() == 3) {
                         PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
                         ghosts.push_back(pjet);
 			}
 		}
 	}
 	
 	Particles subevent = ConstSubPartEvent(origevent, ghosts);
 
 	PseudoJets newevent;
 	foreach (Particle part, subevent) {
 		newevent.push_back(PseudoJet(part.px(), part.py(), part.pz(), part.E()));
 		if (isCharged(part) && part.pT() > _trackptmin) tracks.push_back(part);
 	}
 	
 	JetDefinition jet_def(antikt_algorithm, _jetR);
 	ClusterSequence cs(newevent, jet_def);
 	PseudoJets pjets;
 	pjets = cs.inclusive_jets(10.);
 
 	foreach (PseudoJet pj, pjets) {
 		jetssub.push_back(Jet(pj));
 	}
 	
 #elif MODE==3
       const FinalState fs = applyProjection<FinalState>(event, "FS");
       const Jets loosejets = applyProjection<FastJets>(event, "Jets").jetsByPt(_ptminloose);
 
 	Particles origevent;
 	PseudoJets ghosts;
 
 	foreach (Particle part, fs.particles()) {
 		foreach (Jet jet, loosejets) {
 			if (deltaR(part.momentum(),jet.momentum()) < _dRloose) {
 				origevent.push_back(part);
 				break;
 			}
 			}
 	}
 		
 	foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
 		foreach (Jet jet, loosejets) {
 			FourMomentum mom(p->momentum());
 			if (deltaR(mom,jet.momentum()) < _dRloose && p->status() == 3) {
                         PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
                         ghosts.push_back(pjet);
 				break;
 			}
 		}
 	}
 	
 	Particles subevent = ConstSubPartEvent(origevent, ghosts);
 
 	PseudoJets newevent;
 	foreach (Particle part, subevent) {
 		newevent.push_back(PseudoJet(part.px(), part.py(), part.pz(), part.E()));
 		if (isCharged(part) && (part.pT() > _trackptmin)) tracks.push_back(part);
 	}
 	
 	JetDefinition jet_def(antikt_algorithm, _jetR);
 	ClusterSequence cs(newevent, jet_def);
 	PseudoJets pjets;
 	pjets = cs.inclusive_jets(_ptminloose);
 
 	foreach (PseudoJet pj, pjets) {
 		jetssub.push_back(Jet(pj));
 	}
 
 #endif
       
       
       map<Jet, double, jetcomp> jetmap;
       foreach (Jet jet, jetssub) {
 		double ptreco, ptgen;
 		ptgen = jet.momentum().pT();
 		ptreco = ptsmear(ptgen,cent);
 		jetmap[jet] = ptreco;
       }
 
       Jets jets;
       foreach (Jet jet, jetssub) {
         if (fabs(jet.momentum().eta()) > _etamin && fabs(jet.momentum().eta()) < _etamax && 
 	     jetmap[jet] > _ptmin && jetmap[jet] < _ptmax)
 		jets.push_back(jet);
       }
 
        foreach (Jet jet, jets) {
 		_Njets[0]+=weight;
 		double jetpt = jetmap[jet];
              int ptbin(-1);
              for (size_t l = 0; l < _Nptbins; ++l) {
       	        if (jetpt >= _ptedges[l] && jetpt < _ptedges[l+1]) ptbin=l;
              }
              //shift: first histo is inclusive
              if (ptbin >= 0) ptbin+=1;
              _Njets[ptbin]+=weight;
 		foreach (Particle part, tracks){
 			double dR(deltaR(part,jet));
 			if (dR >= _jetR) continue;
 			double z,xi;
 			z = part.momentum().p3().dot(jet.momentum().p3())/sqr(jet.momentum().p3().mod());
 			//z*=jet.momentum().pT()/jetpt;
 			xi = log(1./z);
 			if (cent < 0.) {
 				if (ptbin > 0) {
 				_h_xi_pp[ptbin]->fill(xi, weight);
 				_h_xi_pp[0]->fill(xi, weight);
 				_h_pt_pp[ptbin]->fill(part.pT(), weight);
 				_h_pt_pp[0]->fill(part.pT(), weight);
 				}
 			}
 		      else {
 				if (ptbin > 0) {
 				_h_xi_PbPb[ptbin]->fill(xi, weight);
 				_h_xi_PbPb[0]->fill(xi, weight);
 				_h_pt_PbPb[ptbin]->fill(part.pT(), weight);
 				_h_pt_PbPb[0]->fill(part.pT(), weight);
 				}
 			}				
 			_h_xi_ratio[ptbin]->fill(xi, weight);
 			_h_xi_ratio[0]->fill(xi, weight);
 			_h_pt_diff[ptbin]->fill(part.pT(), weight);
 			_h_pt_diff[0]->fill(part.pT(), weight);
 		}
 	}
      
 
     }
 
 
     /// Normalise histograms etc., after the run
     void finalize() {
 
 	 for (size_t i = 0; i < _Nptbins+1; ++i) {
 		if (_Njets[i]==0.) _Njets[i]=1.;
 		scale(_h_xi_pp[i],1./_Njets[i]);
 		scale(_h_xi_PbPb[i],1./_Njets[i]);
 		scale(_h_xi_ratio[i],1./_Njets[i]);
 		scale(_h_pt_pp[i],1./_Njets[i]);
 		scale(_h_pt_PbPb[i],1./_Njets[i]);
 		scale(_h_pt_diff[i],1./_Njets[i]);
 	 }
 
     }
 
 
 
   private:
 
     vector<double> _Njets, _ptedges;
     double _jetR, _trackptmin, _etamin, _etamax, _ptmin, _ptmax;
     double _etaloose, _ptminloose, _dRloose;
     size_t _netabins, _nphibins;
     size_t _Nptbins;
 
 
     Histo1DPtr _h_xi_pp[4], _h_xi_PbPb[4], _h_xi_ratio[4];
     Histo1DPtr _h_pt_pp[4], _h_pt_PbPb[4], _h_pt_diff[4];
 
 
   };
 
 
 
   // The hook for the plugin system
 /*#if MODE==0
   DECLARE_RIVET_PLUGIN(CMS_2014_I1299142_4MOMSUB);
 #elif MODE==1
   DECLARE_RIVET_PLUGIN(CMS_2014_I1299142_GRIDSUB2);
 #elif MODE==3
   DECLARE_RIVET_PLUGIN(CMS_2014_I1299142_CONSTSUB);
 #endif*/
   DECLARE_RIVET_PLUGIN(CMS_2014_I1299142);
   
 }
Index: trunk/code/ATLAS_2017_I1511869.cc
===================================================================
--- trunk/code/ATLAS_2017_I1511869.cc	(revision 475)
+++ trunk/code/ATLAS_2017_I1511869.cc	(revision 476)
@@ -1,461 +1,466 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/ChargedFinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
 
 namespace Rivet {
 
 using namespace fastjet;
 
     struct FourMomComp {
         bool operator() (const FourMomentum& p1, const FourMomentum& p2) const
        {return p1.pT()<p2.pT();}
     };
 
 
     struct My4Mom {
 	    double pt;
 	    double mdelta;
 	    double phi;
 	    double y;
     } ;
 
     struct MyPart {
 	    double pt;
 	    double mdelta;
 	    double phi;
 	    double y;
 	    PdgId  id;
     } ;
 
     struct MyDist {
 	    double dR;
 	    size_t ipart;
 	    size_t ighost;
     } ;
     
     bool MyDistComp (const MyDist& dist1, const MyDist& dist2) {
 	    return dist1.dR < dist2.dR; 
     }
 
 
   class ATLAS_2017_I1511869 : public Analysis {
   public:
 
     /// Constructor
     ATLAS_2017_I1511869()
       : Analysis("ATLAS_2017_I1511869")
     {
       _jetR=0.4; _trackptmin=1.; _etamin=0.; _etamax=2.1; _etagapmin=0.8; _etagapmax=1.2; _ptmin=100.; _ptmax=398.; 
 	_Nptbins=4; _Netabins=4;
-	_etaloose=2.8; _ptminloose=70.; _dRloose=0.8;
+	_etaloose=2.5; _ptminloose=90.; _dRloose=1.0;
     }
 
 
     /// @name Analysis methods
     //@{
 
     /// Book histograms and initialise projections before the run
     void init() {
 
- 	 _ptedges += 100.0, 126.0, 158.0, 398.0;
+ 	 _ptedges = {100.0, 126.0, 158.0, 398.0};
       _Njets_pt.resize(_Nptbins, 0.0);
       _Njets_eta.resize(_Netabins, 0.0);
 
 	 FinalState fs(-_etaloose, _etaloose, 0.*GeV);
        addProjection(fs, "FS");
        ChargedFinalState cfs(-_etaloose, _etaloose, _trackptmin*GeV);
        addProjection(cfs, "CFS");
 	
        FastJets fj(fs, FastJets::ANTIKT, _jetR);
        addProjection(fj, "Jets");
 
 	 for (size_t i = 0; i < _Netabins; ++i) {
 		 _h_pt_etabins_pp[i]    = bookHisto1D(i+1, 1, 1);
 		 _h_pt_etabins_PbPb[i]  = bookHisto1D(i+1, 1, 2);
 		 _h_z_etabins_pp[i]     = bookHisto1D(i+5, 1, 1);
 		 _h_z_etabins_PbPb[i]   = bookHisto1D(i+5, 1, 2);
 		 _h_pt_etabins_ratio[i] = bookHisto1D(i*4+9, 1, 1);
 		 _h_z_etabins_ratio[i]  = bookHisto1D(i*4+25, 1, 1);
 	 }
 	 for (size_t i = 0; i < _Nptbins; ++i) {
 		 _h_pt_ptbins_ratio[i] = bookHisto1D(i*4+41, 1, 1);
 		 _h_z_ptbins_ratio[i]  = bookHisto1D(i*4+57, 1, 1);
 	 }
 
     }
 
     
     
     /// extract thermal momenta, that should be subtracted, from HepMC event
     vector<FourMomentum> extractThermalMomenta(const Event & event) {
 	 vector<FourMomentum> thermom;
         foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
             FourMomentum mom(p->momentum());
                if (p->status() == 3) {
                    thermom.push_back(mom);
                }
 	 }
 	 return thermom;
     }
     
     
     /// build map from dummy particles to thermal momenta to speed up subtraction
     map<FourMomentum, FourMomentum, FourMomComp> buildThMomMap(const Event & event, const vector<FourMomentum> * thermom) {
 	    map<FourMomentum, FourMomentum, FourMomComp> thmap;
 	    size_t j(0);
            foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
                   FourMomentum mom(p->momentum());
 			if (fuzzyEquals(mom.E(), 1e-6)) {
 				if (j == thermom->size()) {
 					cout<<"Error: number of dummies does not match number of subtractions, will veto Event"<<endl;
 					thmap.clear();
 					return thmap;
 				}
 				FourMomentum sc((*thermom)[j]);
 				thmap[mom]=sc;
 				j++;
 			}
 	    }
 	    return thmap;
     }
    
     
     /// 4-momentum subtraction
     FourMomentum SubtractJetMom(Jet jet, map<FourMomentum, FourMomentum, FourMomComp> * tmmap) {
         if (tmmap->empty()) return jet.momentum();
         FourMomentum sub(0.,0.,0.,0.);
         foreach (Particle part, jet.constituents()) {
             if (fuzzyEquals(part.E(), 1e-6)) {	 
                     map<FourMomentum, FourMomentum, FourMomComp>::iterator mapit;
                     mapit = (*tmmap).find(part.momentum());
                     if (mapit == (*tmmap).end()) cout<<"Error: did not find matching scattering centre in map.\n"<<part.momentum()<<endl;
                     else sub += mapit->second;
 	     }
         }
         FourMomentum jetmom(jet.momentum());
         //cout<<"Original momentum: "<<jetmom<<endl;
         jetmom -= sub;
         //cout<<"subtracted momentum: "<<jetmom<<endl;
         //cout<<"Subtracted pt: "<<subpt<<"    "<<pjet.pt()<<endl;
         return jetmom;
     }
 
     
 
     /// event-wise Constituent subtraction based on Particles
     Particles ConstSubPartEvent(const Particles particles, PseudoJets pjghosts) {
 	    // sort constituents into two vectors: particles and ghosts (subtraction momenta)
 	    Particles subevent;
 	    vector<MyPart> parts;
 	    vector<My4Mom> ghosts;
 	    
 	    foreach (Particle p, particles) {
 			FourMomentum partmom = p.momentum();
 			if (fuzzyEquals(partmom.E(),1e-6)) continue;
 			if (partmom.E() - fabs(partmom.pz()) > 0.) {
 				MyPart part;
 				part.id = p.pid();
 				part.pt = partmom.pT();
 				part.mdelta = sqrt(partmom.mass2() + sqr(partmom.pT())) - partmom.pT();
 				part.phi = partmom.p3().phi();
 				part.y = partmom.rapidity();
 				//if (abs(part.y) > 5.) {
 				//	cout<<"particle rapidity: "<<part.y<<"     "<<partmom.E()-partmom.pz()<<endl;
 				//	cout<<partmom<<"   "<<part.pt<<endl;
 				//}
 				parts.push_back(part);
 			}
 			else {
 				subevent.push_back(p);
 			}
 	    }
 	    
 	    foreach (PseudoJet pj, pjghosts) {
 			FourMomentum submom = Get4Mom(pj);
 			if (submom.E() - fabs(submom.pz()) > 0.) {
 				My4Mom ghost;
 				ghost.pt = submom.pT();
 				ghost.mdelta = sqrt(submom.mass2() + sqr(submom.pT())) - submom.pT();
 				ghost.phi = submom.p3().phi();
 				ghost.y = submom.rapidity();
 				//cout<<"ghost rapidity: "<<ghost.y<<endl;
 				ghosts.push_back(ghost);
 			}
 			else {
 				PdgId ghostid;
 				if (rand()/RAND_MAX < 1./3.) ghostid = PID::PIPLUS;
 				else if (rand()/RAND_MAX < 2./3.) ghostid = PID::PIMINUS;
 				else ghostid = PID::PI0;
 				Particle pg(ghostid,Get4Mom(pj));
 				subevent.push_back(pg);
 			}
 	    }
 	    
         //cout<<"number of particles: "<<parts.size()<<endl;
         //cout<<"number of ghosts: "<<ghosts.size()<<endl;
 
 	  //create list with all particle-ghosts distances and sort it
 	  vector<MyDist> dists;
 	  for (size_t i=0; i < parts.size(); ++i) {
 		  for (size_t j=0; j < ghosts.size(); ++j) {
 			  //cout<<"i, j = "<<i<<" "<<j<<endl;
 			  MyDist dist;
 			  double Deltaphi(abs(parts[i].phi-ghosts[j].phi));
 			  if (Deltaphi > M_PI) Deltaphi = 2.*M_PI - Deltaphi;
 			  dist.dR = sqrt(sqr(Deltaphi) + sqr(parts[i].y-ghosts[j].y));
 			  dist.ipart = i;
 			  dist.ighost = j;
 			  //cout<<"distance: "<<dist.dR<<endl;
 			  dists.push_back(dist);
 		  }
 	  }
 	  //cout<<"number of pairs: "<<dists.size()<<endl;
 	  //dists.sort(MyDistComp);
 	  std::sort(dists.begin(),dists.end(),MyDistComp);
 	  
 	  //go through all particle-ghost pairs and re-distribute momentum and mass
 	  for (vector<MyDist>::iterator liter=dists.begin(); liter != dists.end(); ++liter) {
 		  //cout<<"dealing with dist "<<liter->dR<<endl;
 		  if (liter->dR > _dRloose) break;
 		  size_t pnum = liter->ipart;
 		  size_t gnum = liter->ighost;
 		  double ptp = parts[pnum].pt;
 		  double ptg = ghosts[gnum].pt;
 		  //cout<<"pts: "<<ptp<<" vs "<<ptg<<endl;
 		  if (ptp > ptg) {
 			  parts[pnum].pt -= ptg;
 			  ghosts[gnum].pt = 0.;
 		  }
 		  else {
 			  ghosts[gnum].pt -= ptp;
 			  parts[pnum].pt = 0.;
 		  }
 		  double mdp = parts[pnum].mdelta;
 		  double mdg = ghosts[gnum].mdelta;
 		  //cout<<"masses: "<<mdp<<" vs "<<mdg<<endl;
 		  if (mdp > mdg) {
 			  parts[pnum].mdelta -= mdg;
 			  ghosts[gnum].mdelta = 0.;
 		  }
 		  else {
 			  ghosts[gnum].mdelta -= mdp;
 			  parts[pnum].mdelta = 0.;
 		  }
 		  
 	  }
 	  
 	  //sum up resulting 4-momenta to get subtracted jet momentum
 	  int nparts(0);
 	  for (size_t i=0; i < parts.size(); ++i) {
 		  if (parts[i].pt > 0.) {
 			  FourMomentum mom((parts[i].pt+parts[i].mdelta)*cosh(parts[i].y), 
 						  parts[i].pt*cos(parts[i].phi),
 						  parts[i].pt*sin(parts[i].phi),
 						 (parts[i].pt+parts[i].mdelta)*sinh(parts[i].y));
 			  //cout<<"particle rapidity: "<<mom.rapidity()<<"  "<<parts[i].y<<endl;
 			  Particle outpart(parts[i].id,mom);
 			  subevent.push_back(outpart);
 			  nparts++;
 		  }
 	  }
 	  int nghost(0);
 	  for (size_t i=0; i < ghosts.size(); ++i) {
 		  if (ghosts[i].pt > 0.) {
 			  FourMomentum mom((ghosts[i].pt+ghosts[i].mdelta)*cosh(ghosts[i].y), 
 						  ghosts[i].pt*cos(ghosts[i].phi),
 						  ghosts[i].pt*sin(ghosts[i].phi),
 						 (ghosts[i].pt+ghosts[i].mdelta)*sinh(ghosts[i].y));
 			  //cout<<"ghost rapidity: "<<mom.rapidity()<<"   "<<ghosts[i].y<<endl;
 			  PdgId ghostid;
 			  if (rand()/RAND_MAX < 1./3.) ghostid = PID::PIPLUS;
 			  else if (rand()/RAND_MAX < 2./3.) ghostid = PID::PIMINUS;
 			  else ghostid = PID::PI0;
 			  Particle outpart(ghostid,mom);
 			  subevent.push_back(outpart);
 			  nghost++;	
 		}
 		  //else cout<<ghosts[i].mdelta<<endl;
 	  }
 	  //cout<<nghost<<" / "<<ghosts.size()+1<<" vs. "<<nparts<<" / "<<parts.size()<<endl;
         return subevent;
     }
 
 
 
 /// 4-momentum adapter
     FourMomentum Get4Mom(PseudoJet pjet){
         FourMomentum mom(pjet.E(), pjet.px(), pjet.py(), pjet.pz());
         return mom;
     }
     
     
     
 
     /// Perform the per-event analysis
     void analyze(const Event& event) {
       const double weight = event.weight();
       const double cent = (event.genEvent()->heavy_ion()?event.genEvent()->heavy_ion()->impact_parameter():-1.);
       if (cent > 0.1) vetoEvent;
       Particles tracks;
 	Jets jetssub;
 
 #if MODE==0      
       const ChargedFinalState cfs = applyProjection<ChargedFinalState>(event, "CFS");
 	tracks = cfs.particles();
       const Jets alljets = applyProjection<FastJets>(event, "Jets").jetsByPt(_ptminloose);
 
       vector<FourMomentum> thermom = extractThermalMomenta(event);
       map<FourMomentum, FourMomentum, FourMomComp> thmommap = buildThMomMap(event, &thermom);
       foreach (Jet jet, alljets) {
 	      FourMomentum jetsub = SubtractJetMom(jet, &thmommap);
 	      jetssub.push_back(Jet(jetsub));
       }
 
 #elif MODE==1
       const FinalState fs = applyProjection<FinalState>(event, "FS");
       const Jets loosejets = applyProjection<FastJets>(event, "Jets").jetsByPt(_ptminloose);
 
 	Particles origevent;
 	PseudoJets ghosts;
-
-	foreach (Jet jet, loosejets) {
 		
-		foreach (Particle part, fs.particles()) {
-			if (deltaR(part.momentum(),jet.momentum()) < _dRloose) origevent.push_back(part);
+	foreach (Particle part, fs.particles()) {
+		foreach (Jet jet, loosejets) {
+			if (deltaR(part.momentum(),jet.momentum()) < _dRloose) {
+				origevent.push_back(part);
+				break;
+			}
 		}
+	}
 		
-		foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
+	foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
+		foreach (Jet jet, loosejets) {
 			FourMomentum mom(p->momentum());
 			if (deltaR(mom,jet.momentum()) < _dRloose && p->status() == 3) {
                         PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
                         ghosts.push_back(pjet);
+				break;
 			}
 		}
 	}
 	
 	Particles subevent = ConstSubPartEvent(origevent, ghosts);
 
 	PseudoJets newevent;
 	foreach (Particle part, subevent) {
 		newevent.push_back(PseudoJet(part.px(), part.py(), part.pz(), part.E()));
 		if (isCharged(part) && (part.pT() > _trackptmin)) tracks.push_back(part);
 	}
 	
 	JetDefinition jet_def(antikt_algorithm, _jetR);
 	ClusterSequence cs(newevent, jet_def);
 	PseudoJets pjets;
 	pjets = cs.inclusive_jets(_ptminloose);
 
 	foreach (PseudoJet pj, pjets) {
 		jetssub.push_back(Jet(pj));
 	}
 
 #endif
 	
        foreach (Jet jet, jetssub) {
 		double jetpt = jet.pT();
 		double jeteta = abs(jet.eta());
 		
 		if (jetpt < _ptmin || jetpt > _ptmax || 
 		    jeteta < _etamin || jeteta > _etamax || 
 		    (jeteta > _etagapmin && jeteta < _etagapmax)) continue;
 
 		int ptbin(-1), etabin(-1);
 		if (jeteta < 0.3) etabin = 1;
 		else if (jeteta < 0.8) etabin = 2;
 		else etabin = 3;
 		
 		for (size_t l = 0; l < _Nptbins; ++l) {
       	        if (jetpt >= _ptedges[l] && jetpt < _ptedges[l+1]) ptbin=l;
 		}
 		//shift: first histo is inclusive
 		if (ptbin >= 0) ptbin+=1;
 		if (ptbin < 0 || etabin < 0) continue;
 		_Njets_eta[0]+=weight;
 		_Njets_eta[etabin]+=weight;
 		_Njets_pt[0]+=weight;
 		_Njets_pt[ptbin]+=weight;
 		
 		foreach (Particle part, tracks){
 			double dR(deltaR(part,jet));
 			if (dR >= _jetR) continue;
 			double z;
 			z = part.pT()*cos(dR)/jetpt;
 
 			if (cent < 0.) {
 				_h_pt_etabins_pp[0]->fill(part.pT(), weight);
 				_h_pt_etabins_pp[etabin]->fill(part.pT(), weight);
 				_h_z_etabins_pp[0]->fill(z, weight);
 				_h_z_etabins_pp[etabin]->fill(z, weight);
 			}
 		      else {
 				_h_pt_etabins_PbPb[0]->fill(part.pT(), weight);
 				_h_pt_etabins_PbPb[etabin]->fill(part.pT(), weight);
 				_h_z_etabins_PbPb[0]->fill(z, weight);
 				_h_z_etabins_PbPb[etabin]->fill(z, weight);
 			}				
 			_h_pt_etabins_ratio[0]->fill(part.pT(), weight);
 			_h_pt_etabins_ratio[etabin]->fill(part.pT(), weight);
 			_h_z_etabins_ratio[0]->fill(z, weight);
 			_h_z_etabins_ratio[etabin]->fill(z, weight);
 			_h_pt_ptbins_ratio[0]->fill(part.pT(), weight);
 			_h_pt_ptbins_ratio[ptbin]->fill(part.pT(), weight);
 			_h_z_ptbins_ratio[0]->fill(z, weight);
 			_h_z_ptbins_ratio[ptbin]->fill(z, weight);
 		}
 	}
      
 
 
     }
 
 
     /// Normalise histograms etc., after the run
     void finalize() {
 
 	 for (size_t i = 0; i < _Netabins; ++i) {
 		if (_Njets_eta[i]==0.) _Njets_eta[i]=1.;
 		scale(_h_pt_etabins_pp[i],1./_Njets_eta[i]);
 		scale(_h_pt_etabins_PbPb[i],1./_Njets_eta[i]);
 		scale(_h_z_etabins_pp[i],1./_Njets_eta[i]);
 		scale(_h_z_etabins_PbPb[i],1./_Njets_eta[i]);
 		scale(_h_pt_etabins_ratio[i],1./_Njets_eta[i]);
 		scale(_h_z_etabins_ratio[i],1./_Njets_eta[i]);
 	 }
 	 for (size_t i = 0; i < _Nptbins; ++i) {
 		if (_Njets_pt[i]==0.) _Njets_pt[i]=1.;
 		scale(_h_pt_ptbins_ratio[i],1./_Njets_pt[i]);
 		scale(_h_z_ptbins_ratio[i],1./_Njets_pt[i]);
 	 }
     }
 
     //@}
 
 
   private:
 
     vector<double> _Njets_pt, _Njets_eta, _ptedges;
     double _jetR, _trackptmin, _etamin, _etamax, _etagapmin, _etagapmax, _ptmin, _ptmax;
     double _etaloose, _ptminloose, _dRloose;
     size_t _Nptbins, _Netabins;
 
 
     Histo1DPtr _h_pt_etabins_pp[4], _h_pt_etabins_PbPb[4];
     Histo1DPtr _h_z_etabins_pp[4], _h_z_etabins_PbPb[4];
     Histo1DPtr _h_pt_etabins_ratio[4];
     Histo1DPtr _h_z_etabins_ratio[4];
     Histo1DPtr _h_pt_ptbins_ratio[4];
     Histo1DPtr _h_z_ptbins_ratio[4];
     
 
   };
 
 
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(ATLAS_2017_I1511869);
 
 
 }