Index: trunk/code/ALICE_2017_I1512107_CONSTSUB.plot
===================================================================
--- trunk/code/ALICE_2017_I1512107_CONSTSUB.plot	(revision 0)
+++ trunk/code/ALICE_2017_I1512107_CONSTSUB.plot	(revision 465)
@@ -0,0 +1,40 @@
+# BEGIN PLOT /ALICE_2017_I1512107_CONSTSUB/chfrac_rat
+LogY=0
+# END PLOT
+
+# BEGIN PLOT /ALICE_2017_I1512107_CONSTSUB/JetMass
+LogY=0
+NormalizeToIntegral=1
+# END PLOT
+
+# BEGIN PLOT /ALICE_2017_I1512107_CONSTSUB/d01-x01-y01
+XLabel=$M_\text{ch jet}\ [\text{GeV}^2]$
+YLabel=$1/N_\text{jets}\ \mathrm{d}N/\mathrm{d}M_\text{ch jet}\ [\mathrm{GeV}^{-2}]$
+Title=$60\,\mathrm{GeV} < p_{\perp,\,\text{ch jet}} < 80\,\mathrm{GeV}$ - 0-10 Pb+Pb
+LogY=0
+NormalizeToIntegral=1
+# END PLOT
+
+# BEGIN PLOT /ALICE_2017_I1512107_CONSTSUB/d02-x01-y01
+XLabel=$M_\text{ch jet}\ [\text{GeV}^2]$
+YLabel=$1/N_\text{jets}\ \mathrm{d}N/\mathrm{d}M_\text{ch jet}\ [\mathrm{GeV}^{-2}]$
+Title=$80\,\mathrm{GeV} < p_{\perp,\,\text{ch jet}} < 100\,\mathrm{GeV}$ - 0-10 Pb+Pb
+LogY=0
+NormalizeToIntegral=1
+# END PLOT
+
+# BEGIN PLOT /ALICE_2017_I1512107_CONSTSUB/d03-x01-y01
+XLabel=$M_\text{ch jet}\ [\text{GeV}^2]$
+YLabel=$1/N_\text{jets}\ \mathrm{d}N/\mathrm{d}M_\text{ch jet}\ [\mathrm{GeV}^{-2}]$
+Title=$100\,\mathrm{GeV} < p_{\perp,\,\text{ch jet}} < 120\,\mathrm{GeV}$ - 0-10 Pb+Pb
+LogY=0
+NormalizeToIntegral=1
+# END PLOT
+
+# BEGIN PLOT /ALICE_2017_I1512107_CONSTSUB/JetpT
+RatioPlotYMin=0
+# END PLOT
+
+# BEGIN PLOT /ALICE_2017_I1512107_CONSTSUB/jetpT
+RatioPlotYMin=0
+# END PLOT
Index: trunk/code/ALICE_CHPART.plot
===================================================================
--- trunk/code/ALICE_CHPART.plot	(revision 464)
+++ trunk/code/ALICE_CHPART.plot	(revision 465)
@@ -1,13 +1,14 @@
 # BEGIN PLOT /ALICE_CHPART/d01-x01-y01
 #Title=[Uncomment and insert title for histogram d01-x01-y01 here]
 #XLabel=[Uncomment and insert x-axis label for histogram d01-x01-y01 here]
 #YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here]
 # + any additional plot settings you might like, see make-plots documentation
 # END PLOT
 
 # ... add more histograms as you need them ...
 
 # BEGIN PLOT /ALICE_CHPART/*
+XMin=2.
 FullRange=1
 # END PLOT
 
Index: trunk/code/ALICE_2017_I1512107.cc
===================================================================
--- trunk/code/ALICE_2017_I1512107.cc	(revision 464)
+++ trunk/code/ALICE_2017_I1512107.cc	(revision 465)
@@ -1,520 +1,719 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/ChargedFinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
 #include "fastjet/ClusterSequence.hh"
 #include "fastjet/ClusterSequence.hh"
 
 #ifndef MODE 
 #warning MODE not defined, using default
 #define MODE 0
 #endif
 
 namespace Rivet {
 
   using namespace fastjet;
+
+    struct PJetComp {
+        bool operator() (const PseudoJet& pj1, const PseudoJet& pj2) const
+       {return pj1.pt()<pj2.pt();}
+    };
 	  
+    struct My4Mom {
+	    double pt;
+	    double mdelta;
+	    double phi;
+	    double y;
+    } ;
+
+    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 ALICE_2017_I1512107_4MOMSUB : public Analysis {
 #elif MODE==1
   class ALICE_2017_I1512107_GRIDSUB1 : public Analysis {
 #elif MODE==2
   class ALICE_2017_I1512107_GRIDSUB2 : public Analysis {
 #elif MODE==3
   class ALICE_2017_I1512107_CHARGED : public Analysis {
+#elif MODE==4
+  class ALICE_2017_I1512107_QSUB : public Analysis {
+#elif MODE==5
+  class ALICE_2017_I1512107_CONSTSUB : public Analysis {
 #endif	  
 	  
   public:
 
     /// Constructor
 #if MODE==0
     ALICE_2017_I1512107_4MOMSUB() : Analysis("ALICE_2017_I1512107_4MOMSUB")
 #elif MODE==1
     ALICE_2017_I1512107_GRIDSUB1() : Analysis("ALICE_2017_I1512107_GRIDSUB1")
 #elif MODE==2
     ALICE_2017_I1512107_GRIDSUB2() : Analysis("ALICE_2017_I1512107_GRIDSUB2")
 #elif MODE==3
     ALICE_2017_I1512107_CHARGED() : Analysis("ALICE_2017_I1512107_CHARGED")
+#elif MODE==4
+    ALICE_2017_I1512107_QSUB() : Analysis("ALICE_2017_I1512107_QSUB")
+#elif MODE==5
+    ALICE_2017_I1512107_CONSTSUB() : Analysis("ALICE_2017_I1512107_CONSTSUB")
 #endif    
     {    
-	    _jetR = 0.4; _ptmin = 50.; _etamax = 0.5; _etaloose = 0.9;
+	    _jetR = 0.4; _ptmin = 20.; _etamax = 0.5; _etaloose = 0.9;
 	    _nphibins = 120; _netabins = 160;
            //_nphibins = 60; _netabins = 80;
 	    _ptscalefac = 0.75; _massscalefac = 0.67;
     }
 
 
     /// @name Analysis methods
     //@{
 
     /// Book histograms and initialise projections before the run
     void init() {
 	    
 	    _ptedges += 60.0, 80.0, 100.0, 120.0;
 
 	    //FinalState fs(-_etaloose, _etaloose, 0.*GeV);
 	    FinalState fs(-3., 3., 0.*GeV);
 	    addProjection(fs,"FS");
 	    ChargedFinalState cfs(-_etaloose, _etaloose, 0.*GeV);
 	    addProjection(cfs,"CFS");
 	    FastJets tfj(cfs, FastJets::ANTIKT, _jetR);
 	    addProjection(tfj, "TrackJets");
 	    FastJets ffj(fs, FastJets::ANTIKT, _jetR);
 	    addProjection(ffj, "FullJets");
 
 
 	    //_h_scalefac = bookHisto1D("scalefac", 50, 0., 1.);
 	    _h_jetmass[0] = bookHisto1D(1, 1, 1);
 	    _h_jetmass[1] = bookHisto1D(2, 1, 1);
 	    _h_jetmass[2] = bookHisto1D(3, 1, 1);
 	    _h_jetpt = bookHisto1D("jetpT", 40, 100., 400.);
 	    _h_jetpt_resc = bookHisto1D("jetpT_resc", 40, 100., 400.);
 	    _h_JetMass = bookHisto1D("JetMass", 50, 0., 100.);
 	    _h_JetpT = bookHisto1D("JetpT", 30, 100, 400);
 	    
 
     }
 
     
     PseudoJets extractThermalMomenta(const Event & event) {
 	 PseudoJets thermom;
         foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
             FourMomentum mom(p->momentum());
                if (p->status() == 3) {
                    PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
                    thermom.push_back(pjet);
                }
 	 }
 	 return thermom;
     }
     
 
+    map<PseudoJet, PseudoJet, PJetComp> buildRecoilMap(const Event & event, const PseudoJets & thermom) {
+        map<PseudoJet, PseudoJet, PJetComp> recmap;
+	  size_t j(0);
+        foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
+            FourMomentum mom(p->momentum());
+            if (mom.E() == 1e-6) {
+	        PseudoJet pj(mom.px(), mom.py(), mom.pz(), mom.E());
+		  if (j == thermom.size()) {
+			cout<<"Error: number of dummies does not match number of subtractions, will veto Event"<<endl;
+			recmap.clear();
+			return recmap;
+		  }
+              PseudoJet sc(thermom[j]);
+              //FourMomentum test(mom);
+              //test *= 10000.;
+	        //FourMomentum test2(Get4Mom(sc));
+		  //double dR(deltaR(test,test2));
+		  //if (dR > 1e-4) continue;
+              recmap[pj]=sc;
+		  j++;
+		}
+	  }
+	  return recmap;
+    }
+
+
     /// 4-momentum subtraction
     FourMomentum SubtractJetMom(Jet jet, PseudoJets * thermom, double fac = 1.) {
         if (thermom->empty()) return jet.momentum();
         FourMomentum sub(0.,0.,0.,0.);
         foreach (Particle part, jet.constituents()) {
             if (part.E() < 1e-5 && part.E() > 1e-7) {
                     FourMomentum test(part.momentum());
                     test *= 10000.;
                     //cout<<"Found dummy :"<<Get4Mom(pj)<<endl;
                     double dRmin(10.);
                     PseudoJet matchthmom;
                     foreach (PseudoJet tm, *thermom) {
                         double dR = deltaR(test.eta(), test.phi(ZERO_2PI), Get4Mom(tm).eta(),Get4Mom(tm).phi(ZERO_2PI));
                         if (dR < dRmin) {
                             dRmin = dR;
                             matchthmom = tm;
                         }
                     }
                     if (dRmin < 1e-5) {
                        //cout<<"Found matching scattering centre :"<<Get4Mom(matchscatcen)<<endl;
                         sub += fac*Get4Mom(matchthmom);
                    }
-                    else cout<<"Error: did not find scattering centre matching dummy.\n";
+                    else cout<<"Error: did not find dummy in map.\n";
             }
         }
         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;
     }
 
+    /// alternative 4-momentum subtraction
+    FourMomentum SubtractJetMom2(PseudoJet pjet, map<PseudoJet, PseudoJet, PJetComp> & recmap) {
+        FourMomentum jetmom(Get4Mom(pjet));
+        if (recmap.size() == 0) return jetmom;
+        FourMomentum sub(0.,0.,0.,0.);
+        foreach (PseudoJet pj, pjet.constituents()) {
+            if (pj.E() == 1e-6) {
+            	map<PseudoJet, PseudoJet, PJetComp>::iterator mapit;
+            	mapit = recmap.find(pj);
+            	if (mapit == recmap.end()) cout<<"Error: did not find matching scattering centre in map.\n";
+			else {
+				FourMomentum submom(Get4Mom(mapit->second));
+				if (submom.mass2() < 0.) sub -= submom;
+				/*{
+					if ((submom.E() < 0. && submom.eta() == -Get4Mom(pj).eta()) || (submom.E() >= 0. && submom.eta() == Get4Mom(pj).eta())) sub -= submom;
+					else sub += submom;
+				}*/
+				else sub += submom;
+			  }
+		} 
+	  }
+        jetmom -= sub;
+        return jetmom;
+    }
+
+
+    /// Constituent subtraction
+    FourMomentum ConstSubJetMom(PseudoJet pjet, map<PseudoJet, PseudoJet, PJetComp> & scmap) {
+	    // sort constituents into two vectors: particles and ghosts (subtraction momenta)
+	    vector<My4Mom> parts;
+	    vector<My4Mom> ghosts;
+	    foreach (PseudoJet pj, pjet.constituents()) {
+            if (fuzzyEquals(pj.E(),1e-6,1e-4)) {
+		      map<PseudoJet, PseudoJet, PJetComp>::iterator mapit;
+			mapit = scmap.find(pj);
+			if (mapit == scmap.end()) cout<<"Error: did not find matching scattering centre in map.\n";
+			else {
+				FourMomentum submom = Get4Mom(mapit->second);
+				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();
+				ghosts.push_back(ghost);
+			}
+		}
+		else {
+			FourMomentum partmom = Get4Mom(pj);
+			My4Mom part;
+			part.pt = partmom.pT();
+			part.mdelta = sqrt(partmom.mass2() + sqr(partmom.pT())) - partmom.pT();
+			part.phi = partmom.p3().phi();
+			part.y = partmom.rapidity();
+			parts.push_back(part);
+		}
+        }
+        //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;
+			  dist.dR = sqrt(sqr(parts[i].phi-ghosts[j].phi) + 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;
+	  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;
+		  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
+	  FourMomentum partmom, ghostmom, jetmom;
+	  for (size_t i=0; i < parts.size(); ++i) {
+		  if (parts[i].pt > 0.) partmom += FourMomentum((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));
+	  }
+	  for (size_t i=0; i < ghosts.size(); ++i) {
+		  if (ghosts[i].pt > 0.) ghostmom += FourMomentum((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<<"particle component: "<<partmom<<endl;
+	  //cout<<"ghost component: "<<ghostmom<<endl;
+	  jetmom = partmom + ghostmom;
+	  //cout<<"jet momentum: "<<jetmom<<endl;
+        return jetmom;
+    }
+
+    
     
 /// initialise and fill grid for grid subtration version 2   
     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. && mom.E() > 1e-5) {
                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] -= mom;
                }
                else {
 		     grid[phibin][etabin] += mom;
                }*/
                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 {
                    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;
     }  
 
     
     /// determine fraction of pt carried by charged particles in jet
     double getScaleFac(const Jet jet, double ptcut, const FinalState * fs = NULL) {
 	    Particles parts;
 	    if (fs) parts = fs->particles();
 	    else parts = jet.constituents();
 	    double scalefac;
 	    double ntracks(0.), nparts(0.);
            FourMomentum chmom(0.,0.,0.,0.), allmom(0.,0.,0.,0.);
 	      foreach (Particle part, parts) {
 		      if (deltaR(part,jet) < _jetR && part.momentum().E() > 1e-5) {
 			      nparts += 1.;
 			      allmom += part.momentum();
                                   if (PID::charge(part.pdgId()) != 0 && part.pT() > ptcut) {
 			                    ntracks += 1.;
 			                    chmom += part.momentum();				      
 			             }
 		      }
 	      }
 	      scalefac = (nparts>0.?ntracks/nparts:1);
 	      //scalefac = (allmom.pT()>0.?chmom.pT()/allmom.pT():1.);
 	    return scalefac;
     }
 
 
 
     /// initialise and fill grid for grid subtration version 1
     void fillJetGrid(Jet jet,  PseudoJets * thermom, vector<vector<FourMomentum> > & grid, vector<vector<FourMomentum> > & thgrid, std::set<std::pair<size_t,size_t> > & gridevent) {
       // fill grid
         foreach (Particle part, jet.constituents()) {
             FourMomentum mom(part.momentum());
                size_t phibin, etabin;
                if (mom.E() < 1e-5 && mom.E() > 1e-7) {
                     /*FourMomentum test(mom);
                     test *= 10000.;
                     //cout<<"Found dummy :"<<Get4Mom(mom)<<endl;
                     double dRmin(10.);
                     PseudoJet matchthmom;
                     foreach (PseudoJet tm, *thermom) {
                         double dR = deltaR(test.eta(), test.phi(ZERO_2PI), Get4Mom(tm).eta(),Get4Mom(tm).phi(ZERO_2PI));
                         if (dR < dRmin) {
                             dRmin = dR;
                             matchthmom = tm;
                         }
                     }
                     if (dRmin < 1e-5) {
 			       mom = Get4Mom(matchthmom);
 				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);}
 				thgrid[phibin][etabin] += mom;*/
 				/*double cellphi((phibin+0.5)*2*M_PI/_nphibins), celleta((etabin+0.5)*8./_netabins-4.);
 				double celltheta(2.*atan(exp(-celleta)));
 				double rho(thgrid[phibin][etabin].E() + mom.E());
 				thgrid[phibin][etabin].setE(rho);
 				thgrid[phibin][etabin].setPx(rho*sin(celltheta)*cos(cellphi));
 				thgrid[phibin][etabin].setPy(rho*sin(celltheta)*sin(cellphi));
 				thgrid[phibin][etabin].setPz(rho*cos(celltheta));*/
 			/*	gridevent.insert(std::pair<size_t,size_t>(phibin,etabin));
                      }
                      else cout<<"Error: did not find scattering centre matching dummy.\n";*/
                }
                else {
 			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);}
 			//grid[phibin][etabin] += mom;
 			double cellphi((phibin+0.5)*2*M_PI/_nphibins), celleta((etabin+0.5)*8./_netabins-4.);
 			double celltheta(2.*atan(exp(-celleta)));
                      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));
 		 }
 	 }
 	 
 	 foreach (PseudoJet pjet, *thermom) {
 		 FourMomentum mom = Get4Mom(pjet);
 		 double dR(deltaR(mom, jet));
 		 if (dR < _jetR) {
 			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);}
 			//thgrid[phibin][etabin] += mom;
 			double cellphi((phibin+0.5)*2*M_PI/_nphibins), celleta((etabin+0.5)*8./_netabins-4.);
 			double celltheta(2.*atan(exp(-celleta)));
                      double rho(thgrid[phibin][etabin].E() + mom.E());
                      thgrid[phibin][etabin].setE(rho);
                      thgrid[phibin][etabin].setPx(rho*sin(celltheta)*cos(cellphi));
                      thgrid[phibin][etabin].setPy(rho*sin(celltheta)*sin(cellphi));
                      thgrid[phibin][etabin].setPz(rho*cos(celltheta));
 			gridevent.insert(std::pair<size_t,size_t>(phibin,etabin));
 		 }
 	 }
 	 
 	 return;
 	 
     }  
 
 
 
 /// 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;
       
-#if MODE==0
+#if (MODE==0 || MODE==4 || MODE==5)
 //       cout<<"Karo:   event number "<<event.genEvent()->event_number()<<endl;
       //Cut cuts = Cuts::abseta < _etaloose && Cuts::pT > _ptmin*GeV;
       Cut cuts = Cuts::abseta < 2.0 && Cuts::pT > _ptmin*GeV;
       const Jets jets = applyProjection<FastJets>(event, "FullJets").jetsByPt(cuts);
 
       PseudoJets thermom = extractThermalMomenta(event);
-      foreach (Jet jet, jets) {
+      map<PseudoJet, PseudoJet, PJetComp> scmap = buildRecoilMap(event, thermom);
+	if (scmap.size() != thermom.size()) {
+		cout<<"Error: size of map does not match number of momenta: "<<scmap.size()<<" vs. "<<thermom.size()<<endl;
+		cout<<"Will veto event."<<endl;
+		vetoEvent;
+	 }
+
+       foreach (Jet jet, jets) {
 	      //double scalefac;
 	      //scalefac = getScaleFac(jet, 0.);
 	      //_h_scalefac->fill(scalefac, weight);
-	      FourMomentum jetsub = SubtractJetMom(jet, &thermom, 1.);
+	      //FourMomentum jetsub = SubtractJetMom(jet, &thermom, 1.);
+#if (MODE==0 || MODE==4)		 
+	      FourMomentum jetsub = SubtractJetMom2(jet, scmap);
+#else
+	      FourMomentum jetsub = ConstSubJetMom(jet, scmap);
+#endif
+		//FourMomentum jetsub = jet.momentum();
 // 	      if (jet.momentum().pT() > 100.) {
 // 		      cout<<"Karo:   "<<jet.momentum().pT()<<"  ->  "<<jetsub.pT()<<endl;
 // 	      }
 	      if (abs(jetsub.eta()) < _etamax) {
 			double mass2(jetsub.mass2());
 			double mass(mass2>0.?sqrt(mass2):-1.);
 			//mass *= scalefac;
 			mass *= _massscalefac;
 			//double jetpt(scalefac*jetsub.pT());
 			double jetpt(jetsub.pT());
 			_h_jetpt->fill(jetpt,weight);
 			jetpt *= _ptscalefac;
 			_h_jetpt_resc->fill(jetpt,weight);			
 			for (size_t i=0; i<3; ++i) {
 				if (jetpt > _ptedges[i] && jetpt <= _ptedges[i+1]) _h_jetmass[i]->fill(mass,weight);
 			}
 	      }
 	      if (jetsub.pT() > 100.*GeV) {
 // 		      cout<<"Karo:   "<<jetsub.pT()<<endl;
 			double mass2(jetsub.mass2());
 			double mass(mass2>0.?sqrt(mass2):-1.);
 			_h_JetMass->fill(mass,weight);
 			_h_JetpT->fill(jetsub.pT(),weight);
 	      }
       }
 #elif MODE==1
       //Cut cuts = Cuts::abseta < _etaloose && Cuts::pT > _ptmin*GeV;
       Cut cuts = Cuts::abseta < 2.0 && Cuts::pT > _ptmin*GeV;
       const Jets jets = applyProjection<FastJets>(event, "FullJets").jetsByPt(cuts);
 
       PseudoJets thermom = extractThermalMomenta(event);
 	    // initialise grid
            vector<vector<FourMomentum> >  grid, thgrid;
            grid.resize(_nphibins);
            thgrid.resize(_nphibins);
            for (size_t i = 0; i < _nphibins; ++i) {
                   grid[i].resize(_netabins);
                   thgrid[i].resize(_netabins);
 		    for (size_t k = 0; k < _nphibins; ++k) {
 			    grid[i][k] = FourMomentum(0., 0., 0., 0.);
 			    thgrid[i][k] = FourMomentum(0., 0., 0., 0.);
 		    }
            }
 	    std::set<std::pair<size_t,size_t> > gridevent;
 
       foreach (Jet jet, jets) {
 		fillJetGrid(jet, &thermom, grid, thgrid, gridevent);
       }
 	 double negenergy(0.);
 	 PseudoJets pevent;
         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;
 	     FourMomentum diff(grid[i][k] - thgrid[i][k]);
             if (grid[i][k].pT() > thgrid[i][k].pT() && diff.E() > 0.) {
                 PseudoJet pjet(diff.px(),diff.py(),diff.pz(),diff.E());
                 pevent.push_back(pjet);
             }
             else negenergy += diff.E();
         }
 	JetDefinition jet_def(antikt_algorithm, _jetR);
 	ClusterSequence cs(pevent, jet_def);
 	PseudoJets pjets = sorted_by_pt(cs.inclusive_jets(_ptmin));
 	//cout<<"size comparison: "<<jets.size()<<"   vs.   "<<pjets.size()<<endl;
       foreach (PseudoJet pjet, pjets) {
 	       FourMomentum jetsub = Get4Mom(pjet);
 	       if (abs(jetsub.eta()) < _etamax) {
 			double mass2(jetsub.mass2());
 			double mass(mass2>0.?sqrt(mass2):-1.);
 			mass *= _massscalefac;
 			double jetpt(jetsub.pT());
 			_h_jetpt->fill(jetpt,weight);
 			jetpt *= _ptscalefac;
 			_h_jetpt_resc->fill(jetpt,weight);			
 			for (size_t i=0; i<3; ++i) {
 				if (jetpt > _ptedges[i] && jetpt <= _ptedges[i+1]){
 					_h_jetmass[i]->fill(mass,weight);
 				}
 			}
 	      }
 	      if (jetsub.pT() > 100.*GeV) {
 			double mass2(jetsub.mass2());
 			double mass(mass2>0.?sqrt(mass2):-1.);
 			_h_JetMass->fill(mass,weight);
 			_h_JetpT->fill(jetsub.pT(),weight);
 	      }
       }
 #elif MODE==2
       const FinalState parts = applyProjection<FinalState>(event, "FS");
 
       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(_ptmin));
       foreach (PseudoJet pjet, pjets) {
 	      FourMomentum jetsub = Get4Mom(pjet);
 	      if (abs(jetsub.eta()) < _etamax) {
 	             //double scalefac;
 		      //scalefac = getScaleFac(Jet(pjet),0.,&parts);
 		      //_h_scalefac->fill(scalefac,weight);
 			double mass2(jetsub.mass2());
 			double mass(mass2>0.?sqrt(mass2):-1.);
 			//mass *= scalefac;
 			mass *= _massscalefac;
 			//double jetpt(scalefac*jetsub.pT());
 			double jetpt(jetsub.pT());
 			_h_jetpt->fill(jetpt,weight);
 			jetpt *= _ptscalefac;
 			_h_jetpt_resc->fill(jetpt,weight);			
 			for (size_t i=0; i<3; ++i) {
 				if (jetpt > _ptedges[i] && jetpt <= _ptedges[i+1]){
 					_h_jetmass[i]->fill(mass,weight);
 				}
 			}
 	      }
 	      if (jetsub.pT() > 100.*GeV) {
 			double mass2(jetsub.mass2());
 			double mass(mass2>0.?sqrt(mass2):-1.);
 			_h_JetMass->fill(mass,weight);
 			_h_JetpT->fill(jetsub.pT(),weight);
 	      }
       }
 #elif MODE==3
       const FinalState parts = applyProjection<FinalState>(event, "FS");
       Cut cuts = Cuts::abseta < _etaloose && Cuts::pT > _ptmin*GeV;
       const Jets jets = applyProjection<FastJets>(event, "TrackJets").jetsByPt(cuts);
 
       foreach (Jet jet, jets) {
 	      //double scalefac;
 	      //scalefac = getScaleFac(jet, 0., &parts);
 	      //_h_scalefac->fill(scalefac, weight);
 	      PseudoJets thermom = extractThermalMomenta(event);
 	      //FourMomentum jetsub = SubtractJetMom(jet, &thermom, scalefac);
 	      FourMomentum jetsub = SubtractJetMom(jet, &thermom, _ptscalefac);
 	      if (abs(jetsub.eta()) < _etamax) {
 			double mass2(jetsub.mass2());
 			double mass(mass2>0.?sqrt(mass2):-1.);
 			double jetpt(jetsub.pT());
 			_h_jetpt_resc->fill(jetpt,weight);			
 			for (size_t i=0; i<3; ++i) {
 				if (jetpt > _ptedges[i] && jetpt <= _ptedges[i+1]) _h_jetmass[i]->fill(mass,weight);
 			}
 	      }
       }
 #endif
 
 
     }
 
 
     /// Normalise histograms etc., after the run
     void finalize() {
 
 	    //normalize(_h_scalefac);
 	    for (size_t i=0; i<3; ++i) {
 		normalize(_h_jetmass[i]);
 	    }
 	    scale(_h_jetpt, crossSection()/sumOfWeights());
 	    scale(_h_jetpt_resc, crossSection()/sumOfWeights());
 	    normalize(_h_JetMass);
 	    scale(_h_JetpT, crossSection()/picobarn/sumOfWeights());
 
     }
 
 
   private:
 
     double _jetR, _ptmin, _etamax, _etaloose;
     double _massscalefac, _ptscalefac;
     size_t _netabins, _nphibins;
     vector<double> _ptedges;
 
 
     //Histo1DPtr _h_scalefac;
     Histo1DPtr _h_jetmass[3];
     Histo1DPtr _h_jetpt, _h_jetpt_resc;
     Histo1DPtr _h_JetMass, _h_JetpT;
     
 
   };
 
 
 
   // The hook for the plugin system
 #if MODE==0
   DECLARE_RIVET_PLUGIN(ALICE_2017_I1512107_4MOMSUB);
 #elif MODE==1
   DECLARE_RIVET_PLUGIN(ALICE_2017_I1512107_GRIDSUB1);
 #elif MODE==2
   DECLARE_RIVET_PLUGIN(ALICE_2017_I1512107_GRIDSUB2);
 #elif MODE==3
   DECLARE_RIVET_PLUGIN(ALICE_2017_I1512107_CHARGED);
+#elif MODE==4
+  DECLARE_RIVET_PLUGIN(ALICE_2017_I1512107_QSUB);
+#elif MODE==5
+  DECLARE_RIVET_PLUGIN(ALICE_2017_I1512107_CONSTSUB);
 #endif
   
 
 
 }
Index: trunk/code/ALICE_2017_I1512107_CONSTSUB.yoda
===================================================================
--- trunk/code/ALICE_2017_I1512107_CONSTSUB.yoda	(revision 0)
+++ trunk/code/ALICE_2017_I1512107_CONSTSUB.yoda	(revision 465)
@@ -0,0 +1,55 @@
+# BEGIN YODA_SCATTER2D /REF/ALICE_2017_I1512107_CONSTSUB/d03-x01-y01
+Path=/REF/ALICE_2017_I1512107_CONSTSUB/d03-x01-y01
+Type=Scatter2D
+# xval	 xerr-	 xerr+	 yval	 yerr-	 yerr+
+1.000000e+00	1.000000e+00	1.000000e+00	2.086430e-03	5.228280e-04	5.228280e-04
+3.000000e+00	1.000000e+00	1.000000e+00	1.294190e-02	1.837800e-03	1.837800e-03
+5.000000e+00	1.000000e+00	1.000000e+00	4.749750e-02	5.917910e-03	5.917910e-03
+7.000000e+00	1.000000e+00	1.000000e+00	7.599270e-02	6.524830e-03	6.524830e-03
+9.000000e+00	1.000000e+00	1.000000e+00	8.285680e-02	6.417800e-03	6.417800e-03
+1.100000e+01	1.000000e+00	1.000000e+00	8.237070e-02	5.847080e-03	5.847080e-03
+1.300000e+01	1.000000e+00	1.000000e+00	6.405460e-02	4.978030e-03	4.978030e-03
+1.500000e+01	1.000000e+00	1.000000e+00	4.619970e-02	3.900820e-03	3.900820e-03
+1.700000e+01	1.000000e+00	1.000000e+00	2.804550e-02	3.062930e-03	3.062930e-03
+1.900000e+01	1.000000e+00	1.000000e+00	2.277020e-02	3.882870e-03	3.882870e-03
+2.100000e+01	1.000000e+00	1.000000e+00	1.399770e-02	3.005680e-03	3.005680e-03
+2.300000e+01	1.000000e+00	1.000000e+00	7.978950e-03	2.001770e-03	2.001770e-03
+2.500000e+01	1.000000e+00	1.000000e+00	9.291080e-03	4.007000e-03	4.007000e-03
+# END YODA_SCATTER2D
+
+# BEGIN YODA_SCATTER2D /REF/ALICE_2017_I1512107_CONSTSUB/d02-x01-y01
+Path=/REF/ALICE_2017_I1512107_CONSTSUB/d02-x01-y01
+Type=Scatter2D
+# xval	 xerr-	 xerr+	 yval	 yerr-	 yerr+
+1.000000e+00	1.000000e+00	1.000000e+00	4.359830e-03	7.934830e-04	7.934830e-04
+3.000000e+00	1.000000e+00	1.000000e+00	3.265260e-02	3.863630e-03	3.863630e-03
+5.000000e+00	1.000000e+00	1.000000e+00	8.170510e-02	6.688520e-03	6.688520e-03
+7.000000e+00	1.000000e+00	1.000000e+00	1.142880e-01	8.765130e-03	8.765130e-03
+9.000000e+00	1.000000e+00	1.000000e+00	9.762640e-02	6.271090e-03	6.271090e-03
+1.100000e+01	1.000000e+00	1.000000e+00	6.621720e-02	3.868830e-03	3.868830e-03
+1.300000e+01	1.000000e+00	1.000000e+00	4.486070e-02	3.064840e-03	3.064840e-03
+1.500000e+01	1.000000e+00	1.000000e+00	2.358060e-02	1.847140e-03	1.847140e-03
+1.700000e+01	1.000000e+00	1.000000e+00	1.451600e-02	1.717860e-03	1.717860e-03
+1.900000e+01	1.000000e+00	1.000000e+00	9.825020e-03	1.861690e-03	1.861690e-03
+2.100000e+01	1.000000e+00	1.000000e+00	7.520780e-03	2.315770e-03	2.315770e-03
+2.300000e+01	1.000000e+00	1.000000e+00	1.568950e-03	4.567960e-04	4.567960e-04
+2.500000e+01	1.000000e+00	1.000000e+00	3.940780e-04	1.718860e-04	1.718860e-04
+# END YODA_SCATTER2D
+
+# BEGIN YODA_SCATTER2D /REF/ALICE_2017_I1512107_CONSTSUB/d01-x01-y01
+Path=/REF/ALICE_2017_I1512107_CONSTSUB/d01-x01-y01
+Type=Scatter2D
+# xval	 xerr-	 xerr+	 yval	 yerr-	 yerr+
+1.000000e+00	1.000000e+00	1.000000e+00	4.902490e-03	4.382700e-04	4.382700e-04
+3.000000e+00	1.000000e+00	1.000000e+00	4.526230e-02	3.077680e-03	3.077680e-03
+5.000000e+00	1.000000e+00	1.000000e+00	1.138460e-01	6.931320e-03	6.931320e-03
+7.000000e+00	1.000000e+00	1.000000e+00	1.342230e-01	5.855050e-03	5.855050e-03
+9.000000e+00	1.000000e+00	1.000000e+00	9.829660e-02	3.887810e-03	3.887810e-03
+1.100000e+01	1.000000e+00	1.000000e+00	6.168980e-02	2.436980e-03	2.436980e-03
+1.300000e+01	1.000000e+00	1.000000e+00	2.564310e-02	1.312140e-03	1.312140e-03
+1.500000e+01	1.000000e+00	1.000000e+00	1.025280e-02	7.268920e-04	7.268920e-04
+1.700000e+01	1.000000e+00	1.000000e+00	4.897520e-03	9.256550e-04	9.256550e-04
+1.900000e+01	1.000000e+00	1.000000e+00	5.489800e-04	9.620230e-05	9.620230e-05
+2.100000e+01	1.000000e+00	1.000000e+00	4.355450e-04	2.104840e-04	2.104840e-04
+# END YODA_SCATTER2D
+
Index: trunk/code/ALICE_2017_I1512107_CONSTSUB.cc
===================================================================
--- trunk/code/ALICE_2017_I1512107_CONSTSUB.cc	(revision 0)
+++ trunk/code/ALICE_2017_I1512107_CONSTSUB.cc	(revision 465)
@@ -0,0 +1,934 @@
+// -*- C++ -*-
+#include "Rivet/Analysis.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Projections/ChargedFinalState.hh"
+#include "Rivet/Projections/FastJets.hh"
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/ClusterSequence.hh"
+
+namespace Rivet {
+
+  using namespace fastjet;
+
+    struct PJetComp {
+        bool operator() (const PseudoJet& pj1, const PseudoJet& pj2) const
+       {return pj1.pt()<pj2.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; 
+    }
+
+
+//#if MODE==0
+  class ALICE_2017_I1512107_CONSTSUB : public Analysis {
+//#endif	  
+	  
+  public:
+
+    /// Constructor
+//#if MODE==0
+    ALICE_2017_I1512107_CONSTSUB() : Analysis("ALICE_2017_I1512107_CONSTSUB")
+//#endif    
+    {    
+	    _jetR = 0.4; _ptmin = 60.; _etamax = 0.5; _etaloose = 1.5;
+	    _ptscalefac = 0.75; _massscalefac = 0.67;
+    }
+
+
+    /// @name Analysis methods
+    //@{
+
+    /// Book histograms and initialise projections before the run
+    void init() {
+	    
+	    _ptedges += 60.0, 80.0, 100.0, 120.0;
+
+	    FinalState fs(-_etaloose, _etaloose, 0.*GeV);
+	    //FinalState fs(-3., 3., 0.*GeV);
+	    addProjection(fs,"FS");
+	    //ChargedFinalState cfs(-_etaloose, _etaloose, 0.*GeV);
+	    //addProjection(cfs,"CFS");
+	    //FastJets tfj(cfs, FastJets::ANTIKT, _jetR);
+	    //addProjection(tfj, "TrackJets");
+	    FastJets ffj(fs, FastJets::ANTIKT, _jetR);
+	    addProjection(ffj, "FullJets");
+
+
+	    //_h_scalefac = bookHisto1D("scalefac", 50, 0., 1.);
+	    _h_jetmass[0] = bookHisto1D(1, 1, 1);
+	    _h_jetmass[1] = bookHisto1D(2, 1, 1);
+	    _h_jetmass[2] = bookHisto1D(3, 1, 1);
+	    _h_jetpt = bookHisto1D("jetpT", 40, 100., 400.);
+	    _h_jetpt_resc = bookHisto1D("jetpT_resc", 40, 100., 400.);
+	    _h_JetMass = bookHisto1D("JetMass", 50, 0., 100.);
+	    _h_JetpT = bookHisto1D("JetpT", 30, 100, 400);
+	    
+
+    }
+
+    
+    
+    /// event-wise Constituent subtraction
+    PseudoJets ConstSubEvent(PseudoJets pjparts, PseudoJets pjghosts) {
+	    // sort constituents into two vectors: particles and ghosts (subtraction momenta)
+	    PseudoJets subevent;
+	    vector<My4Mom> parts;
+	    vector<My4Mom> ghosts;
+	    
+	    foreach (PseudoJet pj, pjparts) {
+			FourMomentum partmom = Get4Mom(pj);
+			if (partmom.E() - fabs(partmom.pz()) > 0.) {
+				My4Mom part;
+				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(pj);
+			}
+	    }
+	    
+	    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 {
+				subevent.push_back(pj);
+			}
+	    }
+	    
+        //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 > 0.3) 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;
+                    PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+			  subevent.push_back(pjet);
+			  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;
+                    PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+			  subevent.push_back(pjet);
+			  nghost++;	
+		}
+		  //else cout<<ghosts[i].mdelta<<endl;
+	  }
+	  cout<<nghost<<" / "<<ghosts.size()+1<<" vs. "<<nparts<<" / "<<parts.size()<<endl;
+        return subevent;
+    }
+
+
+
+    /// 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 {
+				cout<<p.pid()<<"   "<<partmom<<endl;
+				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 {
+				cout<<submom<<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 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 > 0.5) 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;
+    }
+
+
+
+    /// jet-wise Constituent subtraction
+    FourMomentum ConstSubJetMom(Jet jet, map<PseudoJet, PseudoJet, PJetComp> & scmap) {
+	    // sort constituents into two vectors: particles and ghosts (subtraction momenta)
+	    vector<My4Mom> parts;
+	    vector<My4Mom> ghosts;
+	    foreach (Particle p, jet.constituents()) {
+		PseudoJet pj(p.px(), p.py(), p.pz(), p.E());    
+            if (fuzzyEquals(pj.E(),1e-6,1e-4)) {
+		      map<PseudoJet, PseudoJet, PJetComp>::iterator mapit;
+			mapit = scmap.find(pj);
+			if (mapit == scmap.end()) cout<<"Error: did not find matching scattering centre in map.\n";
+			else {
+				FourMomentum submom = Get4Mom(mapit->second);
+				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();
+				ghosts.push_back(ghost);
+			}
+		}
+		else {
+			FourMomentum partmom = Get4Mom(pj);
+			My4Mom part;
+			part.pt = partmom.pT();
+			part.mdelta = sqrt(partmom.mass2() + sqr(partmom.pT())) - partmom.pT();
+			part.phi = partmom.p3().phi();
+			part.y = partmom.rapidity();
+			parts.push_back(part);
+		}
+        }
+        //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;
+	  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;
+		  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
+	  FourMomentum partmom, ghostmom, jetmom;
+	  int nparts(0);
+	  for (size_t i=0; i < parts.size(); ++i) {
+		  if (parts[i].pt > 0.) {
+			  partmom += FourMomentum((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));
+			  nparts++;
+		  }
+	  }
+	  int 	nghost(0);
+	  for (size_t i=0; i < ghosts.size(); ++i) {
+		  if (ghosts[i].pt > 0.) {
+			  ghostmom += FourMomentum((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));
+			  nghost++;
+		  }
+	  }
+	  //cout<<"particle component: "<<partmom<<endl;
+	  //cout<<"ghost component: "<<ghostmom<<endl;
+	  jetmom = partmom + ghostmom;
+	  //cout<<"jet momentum: "<<jetmom<<endl;
+	  if (jetmom.pT() > 60. && nghost > 0) cout<<nghost<<" / "<<ghosts.size()+1<<" vs. "<<nparts<<" / "<<parts.size()<<"  for "<<jetmom.pT()<<endl;
+        return jetmom;
+    }
+
+    /// jet-wise Constituent subtraction
+    Particles ConstSubPartJetMom(Jet jet, map<PseudoJet, PseudoJet, PJetComp> & scmap) {
+	    // sort constituents into two vectors: particles and ghosts (subtraction momenta)
+	    vector<MyPart> parts;
+	    vector<My4Mom> ghosts;
+	    foreach (Particle p, jet.constituents()) {
+            if (fuzzyEquals(p.E(),1e-6,1e-4)) {
+			PseudoJet pj(p.px(), p.py(), p.pz(), p.E());
+		      map<PseudoJet, PseudoJet, PJetComp>::iterator mapit;
+			mapit = scmap.find(pj);
+			if (mapit == scmap.end()) cout<<"Error: did not find matching scattering centre in map.\n";
+			else {
+				FourMomentum submom = Get4Mom(mapit->second);
+				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();
+				ghosts.push_back(ghost);
+			}
+		}
+		else {
+			FourMomentum partmom = p.momentum();
+			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();
+			parts.push_back(part);
+		}
+        }
+        //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;
+	  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;
+		  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.;
+		  }
+	  }
+
+	  Particles subevent;
+	  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);
+		  }
+	  }
+	  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);
+		}
+	  }
+        return subevent;
+    }
+    
+     PseudoJets extractThermalMomenta(const Event & event) {
+	 PseudoJets thermom;
+        foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
+            FourMomentum mom(p->momentum());
+               if (p->status() == 3) {
+                   PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+                   thermom.push_back(pjet);
+               }
+	 }
+	 return thermom;
+    }
+    
+
+    map<PseudoJet, PseudoJet, PJetComp> buildRecoilMap(const Event & event, const PseudoJets & thermom) {
+        map<PseudoJet, PseudoJet, PJetComp> recmap;
+	  size_t j(0);
+        foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
+            FourMomentum mom(p->momentum());
+            if (mom.E() == 1e-6) {
+	        PseudoJet pj(mom.px(), mom.py(), mom.pz(), mom.E());
+		  if (j == thermom.size()) {
+			cout<<"Error: number of dummies does not match number of subtractions, will veto Event"<<endl;
+			recmap.clear();
+			return recmap;
+		  }
+              PseudoJet sc(thermom[j]);
+              //FourMomentum test(mom);
+              //test *= 10000.;
+	        //FourMomentum test2(Get4Mom(sc));
+		  //double dR(deltaR(test,test2));
+		  //if (dR > 1e-4) continue;
+              recmap[pj]=sc;
+		  j++;
+		}
+	  }
+	  return recmap;
+    }
+
+
+   
+/// determine fraction of pt carried by charged particles in jet
+    double getScaleFac(const Jet jet, double ptcut, const FinalState * fs = NULL) {
+	    Particles parts;
+	    if (fs) parts = fs->particles();
+	    else parts = jet.constituents();
+	    double scalefac;
+	    double ntracks(0.), nparts(0.);
+           FourMomentum chmom(0.,0.,0.,0.), allmom(0.,0.,0.,0.);
+	      foreach (Particle part, parts) {
+		      if (deltaR(part,jet) < _jetR && part.momentum().E() > 1e-5) {
+			      nparts += 1.;
+			      allmom += part.momentum();
+                                  if (PID::charge(part.pdgId()) != 0 && part.pT() > ptcut) {
+			                    ntracks += 1.;
+			                    chmom += part.momentum();				      
+			             }
+		      }
+	      }
+	      scalefac = (nparts>0.?ntracks/nparts:1);
+	      //scalefac = (allmom.pT()>0.?chmom.pT()/allmom.pT():1.);
+	    return scalefac;
+    }
+
+
+
+
+
+/// 4-momentum adapter
+    FourMomentum Get4Mom(PseudoJet pjet){
+        FourMomentum mom(pjet.E(), pjet.px(), pjet.py(), pjet.pz());
+        return mom;
+    }
+    
+
+/// PseudoJet adapter
+    PseudoJet GetPseudoJet(FourMomentum mom){
+        PseudoJet pjet(mom.px(), mom.py(), mom.pz(), mom.E());
+        return pjet;
+    }
+    
+    
+    
+
+    /// 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;
+      
+#if (MODE==0)
+// jet-wise subtraction	based on momenta
+      //Cut cuts = Cuts::abseta < _etaloose && Cuts::pT > _ptmin*GeV;
+      Cut cuts = Cuts::abseta < 2.0 && Cuts::pT > _ptmin*GeV;
+      const Jets jets = applyProjection<FastJets>(event, "FullJets").jetsByPt(cuts);
+
+      PseudoJets thermom = extractThermalMomenta(event);
+      map<PseudoJet, PseudoJet, PJetComp> scmap = buildRecoilMap(event, thermom);
+	if (scmap.size() != thermom.size()) {
+		cout<<"Error: size of map does not match number of momenta: "<<scmap.size()<<" vs. "<<thermom.size()<<endl;
+		cout<<"Will veto event."<<endl;
+		vetoEvent;
+	 }
+
+       foreach (Jet jet, jets) {
+	      //double scalefac;
+	      //scalefac = getScaleFac(jet, 0.);
+	      //_h_scalefac->fill(scalefac, weight);
+	      //FourMomentum jetsub = SubtractJetMom(jet, &thermom, 1.);
+	      FourMomentum jetsub = ConstSubJetMom(jet, scmap);
+	      if (abs(jetsub.eta()) < _etamax) {
+			double mass2(jetsub.mass2());
+			double mass(mass2>0.?sqrt(mass2):-1.);
+			//mass *= scalefac;
+			mass *= _massscalefac;
+			//double jetpt(scalefac*jetsub.pT());
+			double jetpt(jetsub.pT());
+			_h_jetpt->fill(jetpt,weight);
+			jetpt *= _ptscalefac;
+			_h_jetpt_resc->fill(jetpt,weight);			
+			for (size_t i=0; i<3; ++i) {
+				if (jetpt > _ptedges[i] && jetpt <= _ptedges[i+1]) _h_jetmass[i]->fill(mass,weight);
+			}
+	      }
+	      if (jetsub.pT() > 100.*GeV) {
+// 		      cout<<"Karo:   "<<jetsub.pT()<<endl;
+			double mass2(jetsub.mass2());
+			double mass(mass2>0.?sqrt(mass2):-1.);
+			_h_JetMass->fill(mass,weight);
+			_h_JetpT->fill(jetsub.pT(),weight);
+	      }
+      }
+      
+#elif (MODE==1)
+// event-wise subtraction based on momenta
+
+	PseudoJets parts, ghosts, pevent;
+
+	foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
+            FourMomentum mom(p->momentum());
+		if (fabs(mom.eta()) < _etaloose) {
+	      if(p->status() == 1) {
+                   PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+			 if (mom.E() > 1e-6) parts.push_back(pjet);
+			 pevent.push_back(pjet);
+		}
+		else if (p->status() == 3) {
+                   PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+                   ghosts.push_back(pjet);
+		}
+		else if(p->status() == 4) {
+                   PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+			 parts.push_back(pjet);
+		}
+		}
+	}
+	  
+	PseudoJets newevent;
+	newevent = ConstSubEvent(parts, ghosts);
+	
+	
+	JetDefinition jet_def(antikt_algorithm, _jetR);
+	ClusterSequence cs(newevent, jet_def);
+	PseudoJets pjets;
+	pjets = sorted_by_pt(cs.inclusive_jets(_ptmin));
+	
+	
+	foreach (PseudoJet pj, pjets) {
+		FourMomentum jet(Get4Mom(pj));
+		
+		if (fabs(jet.eta()) < _etamax) {
+			double mass2(jet.mass2());
+			double mass(mass2>0.?sqrt(mass2):-1.);
+			//mass *= scalefac;
+			mass *= _massscalefac;
+			//double jetpt(scalefac*jet.pT());
+			double jetpt(jet.pT());
+			_h_jetpt->fill(jetpt,weight);
+			jetpt *= _ptscalefac;
+			_h_jetpt_resc->fill(jetpt,weight);			
+			for (size_t i=0; i<3; ++i) {
+				if (jetpt > _ptedges[i] && jetpt <= _ptedges[i+1]) _h_jetmass[i]->fill(mass,weight);
+			}
+	      }
+	      if (jet.pT() > 100.*GeV) {
+// 		      cout<<"Karo:   "<<jetsub.pT()<<endl;
+			double mass2(jet.mass2());
+			double mass(mass2>0.?sqrt(mass2):-1.);
+			_h_JetMass->fill(mass,weight);
+			_h_JetpT->fill(jet.pT(),weight);
+	      }
+      }
+      
+#elif (MODE==2)
+// jet-wise subtraction	based on particles
+
+      Cut cuts = Cuts::abseta < 2.0 && Cuts::pT > _ptmin*GeV;
+      const Jets jets = applyProjection<FastJets>(event, "FullJets").jetsByPt(cuts);
+
+      PseudoJets thermom = extractThermalMomenta(event);
+      map<PseudoJet, PseudoJet, PJetComp> scmap = buildRecoilMap(event, thermom);
+	if (scmap.size() != thermom.size()) {
+		cout<<"Error: size of map does not match number of momenta: "<<scmap.size()<<" vs. "<<thermom.size()<<endl;
+		cout<<"Will veto event."<<endl;
+		vetoEvent;
+	 }
+
+       foreach (Jet jet, jets) {
+	      Particles jetpartsub = ConstSubPartJetMom(jet, scmap);
+		FourMomentum jetsub, fulljet;
+		foreach (Particle p, jetpartsub) {
+			if (isCharged(p)) jetsub += p.momentum();
+			fulljet += p.momentum();
+		}
+	      if (abs(jetsub.eta()) < _etamax) {
+			double mass2(jetsub.mass2());
+			double mass(mass2>0.?sqrt(mass2):-1.);
+			double jetpt(jetsub.pT());
+			_h_jetpt_resc->fill(jetpt,weight);			
+			for (size_t i=0; i<3; ++i) {
+				if (jetpt > _ptedges[i] && jetpt <= _ptedges[i+1]) _h_jetmass[i]->fill(mass,weight);
+			}
+	      }
+
+	      if (abs(fulljet.eta()) < _etamax) {
+			double jetpt(fulljet.pT());
+			_h_jetpt->fill(jetpt,weight);
+	      }
+	      if (fulljet.pT() > 100.*GeV) {
+			double mass2(fulljet.mass2());
+			double mass(mass2>0.?sqrt(mass2):-1.);
+			_h_JetMass->fill(mass,weight);
+			_h_JetpT->fill(fulljet.pT(),weight);
+	      }
+      }
+      
+
+#else 
+// event-wise subtraction based on particles
+
+      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) {
+		if (isCharged(part)) newevent.push_back(PseudoJet(part.px(), part.py(), part.pz(), part.E()));
+	}
+	
+	JetDefinition jet_def(antikt_algorithm, _jetR);
+	ClusterSequence cs(newevent, jet_def);
+	PseudoJets pjets;
+	pjets = sorted_by_pt(cs.inclusive_jets(_ptmin));
+
+	ClusterSequence fullcs(subevent, jet_def);
+	PseudoJets fullpjets;
+	fullpjets = sorted_by_pt(fullcs.inclusive_jets(_ptmin));
+
+		
+	
+	foreach (PseudoJet pj, pjets) {
+		FourMomentum jet(Get4Mom(pj));
+		
+		if (fabs(jet.eta()) < _etamax) {
+			double mass2(jet.mass2());
+			double mass(mass2>0.?sqrt(mass2):-1.);
+			double jetpt(jet.pT());
+			_h_jetpt_resc->fill(jetpt,weight);			
+			for (size_t i=0; i<3; ++i) {
+				if (jetpt > _ptedges[i] && jetpt <= _ptedges[i+1]) _h_jetmass[i]->fill(mass,weight);
+			}
+	      }
+      }
+	
+	foreach (PseudoJet pj, fullpjets) {
+		FourMomentum jet(Get4Mom(pj));
+		
+		if (fabs(jet.eta()) < _etamax) {
+			double jetpt(jet.pT());
+			_h_jetpt->fill(jetpt,weight);
+	      }
+	      if (jet.pT() > 100.*GeV) {
+			double mass2(jet.mass2());
+			double mass(mass2>0.?sqrt(mass2):-1.);
+			_h_JetMass->fill(mass,weight);
+			_h_JetpT->fill(jet.pT(),weight);
+	      }
+      }
+#endif
+      
+
+    }
+
+
+    /// Normalise histograms etc., after the run
+    void finalize() {
+
+	    //normalize(_h_scalefac);
+	    for (size_t i=0; i<3; ++i) {
+		normalize(_h_jetmass[i]);
+	    }
+	    scale(_h_jetpt, crossSection()/sumOfWeights());
+	    scale(_h_jetpt_resc, crossSection()/sumOfWeights());
+	    normalize(_h_JetMass);
+	    scale(_h_JetpT, crossSection()/picobarn/sumOfWeights());
+
+    }
+
+
+  private:
+
+    double _jetR, _ptmin, _etamax, _etaloose;
+    double _massscalefac, _ptscalefac;
+    size_t _netabins, _nphibins;
+    vector<double> _ptedges;
+
+
+    //Histo1DPtr _h_scalefac;
+    Histo1DPtr _h_jetmass[3];
+    Histo1DPtr _h_jetpt, _h_jetpt_resc;
+    Histo1DPtr _h_JetMass, _h_JetpT;
+    
+
+  };
+
+
+
+  // The hook for the plugin system
+//#if MODE==0
+  DECLARE_RIVET_PLUGIN(ALICE_2017_I1512107_CONSTSUB);
+//#endif
+  
+
+
+}
Index: trunk/code/ALICE_2017_I1512107_4MOMSUB.plot
===================================================================
--- trunk/code/ALICE_2017_I1512107_4MOMSUB.plot	(revision 464)
+++ trunk/code/ALICE_2017_I1512107_4MOMSUB.plot	(revision 465)
@@ -1,28 +1,40 @@
 # BEGIN PLOT /ALICE_2017_I1512107_4MOMSUB/chfrac_rat
 LogY=0
 # END PLOT
 
-# BEGIN PLOT /ALICE_2017_I1512107_4MOMSUB/jetmass
+# BEGIN PLOT /ALICE_2017_I1512107_4MOMSUB/JetMass
 LogY=0
+NormalizeToIntegral=1
 # END PLOT
 
 # BEGIN PLOT /ALICE_2017_I1512107_4MOMSUB/d01-x01-y01
 XLabel=$M_\text{ch jet}\ [\text{GeV}^2]$
 YLabel=$1/N_\text{jets}\ \mathrm{d}N/\mathrm{d}M_\text{ch jet}\ [\mathrm{GeV}^{-2}]$
 Title=$60\,\mathrm{GeV} < p_{\perp,\,\text{ch jet}} < 80\,\mathrm{GeV}$ - 0-10 Pb+Pb
 LogY=0
+NormalizeToIntegral=1
 # END PLOT
 
 # BEGIN PLOT /ALICE_2017_I1512107_4MOMSUB/d02-x01-y01
 XLabel=$M_\text{ch jet}\ [\text{GeV}^2]$
 YLabel=$1/N_\text{jets}\ \mathrm{d}N/\mathrm{d}M_\text{ch jet}\ [\mathrm{GeV}^{-2}]$
 Title=$80\,\mathrm{GeV} < p_{\perp,\,\text{ch jet}} < 100\,\mathrm{GeV}$ - 0-10 Pb+Pb
 LogY=0
+NormalizeToIntegral=1
 # END PLOT
 
 # BEGIN PLOT /ALICE_2017_I1512107_4MOMSUB/d03-x01-y01
 XLabel=$M_\text{ch jet}\ [\text{GeV}^2]$
 YLabel=$1/N_\text{jets}\ \mathrm{d}N/\mathrm{d}M_\text{ch jet}\ [\mathrm{GeV}^{-2}]$
 Title=$100\,\mathrm{GeV} < p_{\perp,\,\text{ch jet}} < 120\,\mathrm{GeV}$ - 0-10 Pb+Pb
 LogY=0
+NormalizeToIntegral=1
+# END PLOT
+
+# BEGIN PLOT /ALICE_2017_I1512107_4MOMSUB/JetpT
+RatioPlotYMin=0
+# END PLOT
+
+# BEGIN PLOT /ALICE_2017_I1512107_4MOMSUB/jetpT
+RatioPlotYMin=0
 # END PLOT