SISConePlugin.cc

Go to the documentation of this file.
00001 
00002 // fastjet stuff
00003 #include "fastjet/ClusterSequence.hh"
00004 #include "fastjet/SISConePlugin.hh"
00005 
00006 // scones stuff
00007 #include "siscone/momentum.h"
00008 #include "siscone/siscone.h"
00009 
00010 // other stuff
00011 #include<sstream>
00012 
00013 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00014 
00015 using namespace std;
00016 using namespace siscone;
00017 
00018 // the static objects
00019 auto_ptr<SISConePlugin     > SISConePlugin::stored_plugin    ;
00020 auto_ptr<vector<PseudoJet> > SISConePlugin::stored_particles ;
00021 auto_ptr<Csiscone          > SISConePlugin::stored_siscone   ;
00022 
00023 string SISConePlugin::description () const {
00024   ostringstream desc;
00025   
00026   const string on = "on";
00027   const string off = "off";
00028 
00029   string sm_scale_string = "split-merge uses " + 
00030     split_merge_scale_name(Esplit_merge_scale(split_merge_scale()));
00031 
00032   desc << "SISCone jet algorithm with " ;
00033   desc << "cone_radius = "       << cone_radius        () << ", ";
00034   desc << "overlap_threshold = " << overlap_threshold  () << ", ";
00035   desc << "n_pass_max = "        << n_pass_max         () << ", ";
00036   desc << "protojet_ptmin = "    << protojet_ptmin()      << ", ";
00037   desc <<  sm_scale_string                                << ", ";
00038   desc << "caching turned "      << (caching() ? on : off);
00039   desc << ", SM stop scale = "     << _split_merge_stopping_scale;
00040 
00041 
00042   // create a fake siscone object so that we can find out more about it
00043   Csiscone siscone;
00044   if (siscone.merge_identical_protocones) {
00045     desc << ", and (IR unsafe) merge_indentical_protocones=true" ;
00046   }
00047 
00048   desc << ", SISCone code v" << VERSION;
00049 
00050   return desc.str();
00051 }
00052 
00053 
00055 template<> PseudoJet::PseudoJet(const Cmomentum & four_vector) {
00056   (*this) = PseudoJet(four_vector.px,four_vector.py,four_vector.pz,
00057                       four_vector.E);
00058 }
00059 
00060 
00061 void SISConePlugin::run_clustering(ClusterSequence & clust_seq) const {
00062 
00063 
00064   Csiscone * siscone;
00065   Csiscone   local_siscone;
00066 
00067   unsigned n = clust_seq.jets().size();
00068 
00069   bool new_siscone = true; // by default we'll be running it
00070 
00071   if (caching()) {
00072 
00073     // Establish if we have a cached run with the same R, npass and
00074     // particles. If not then do any tidying up / reallocation that's
00075     // necessary for the next round of caching, otherwise just set
00076     // relevant pointers so that we can reuse and old run.
00077     if (stored_siscone.get() != 0) {
00078       new_siscone = !(stored_plugin->cone_radius()   == cone_radius()
00079                       && stored_plugin->n_pass_max() == n_pass_max()  
00080                       && stored_particles->size()    == n);
00081       if (!new_siscone) {
00082         for(unsigned i = 0; i < n; i++) {
00083           // only check momentum because indices will be correctly dealt
00084           // with anyway when extracting the clustering order.
00085           new_siscone |= !have_same_momentum(clust_seq.jets()[i], 
00086                                              (*stored_particles)[i]);
00087         }
00088       }
00089     } 
00090       
00091     // allocate the new siscone, etc., if need be
00092     if (new_siscone) {
00093       stored_siscone  .reset( new Csiscone                           );
00094       stored_particles.reset( new vector<PseudoJet>(clust_seq.jets()));
00095       stored_plugin   .reset( new SISConePlugin(*this)               );
00096     }
00097 
00098     siscone = stored_siscone.get();
00099   } else {
00100     siscone = &local_siscone;
00101   }
00102 
00103   // make sure stopping scale is set in siscone
00104   siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale;
00105   // when running with ghosts for passive areas, do not put the
00106   // ghosts into the stable-cone search (not relevant)
00107   siscone->stable_cone_soft_pt2_cutoff = ghost_separation_scale()
00108                                          * ghost_separation_scale();
00109 
00110   if (new_siscone) {
00111     // transfer fastjet initial particles into the siscone type
00112     vector<Cmomentum> siscone_momenta(n);
00113     for(unsigned i = 0; i < n; i++) {
00114       const PseudoJet & p = clust_seq.jets()[i]; // shorthand
00115       siscone_momenta[i] = Cmomentum(p.px(), p.py(), p.pz(), p.E());
00116     }
00117     
00118     // run the jet finding
00119     //cout << "plg sms: " << split_merge_scale() << endl;
00120     siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
00121                           n_pass_max(), protojet_or_ghost_ptmin(), 
00122                           Esplit_merge_scale(split_merge_scale()));
00123   } else {
00124     // just run the overlap part of the jets.
00125     //cout << "plg rcmp sms: " << split_merge_scale() << endl;
00126     siscone->recompute_jets(overlap_threshold(), protojet_or_ghost_ptmin(), 
00127                             Esplit_merge_scale(split_merge_scale()));
00128   }
00129 
00130   // extract the jets [in reverse order -- to get nice ordering in pt at end]
00131   int njet = siscone->jets.size();
00132 
00133   for (int ijet = njet-1; ijet >= 0; ijet--) {
00134     const Cjet & jet = siscone->jets[ijet]; // shorthand
00135     
00136     // Successively merge the particles that make up the cone jet
00137     // until we have all particles in it.  Start off with the zeroth
00138     // particle.
00139     int jet_k = jet.contents[0];
00140     for (unsigned ipart = 1; ipart < jet.contents.size(); ipart++) {
00141       // take the last result of the merge
00142       int jet_i = jet_k;
00143       // and the next element of the jet
00144       int jet_j = jet.contents[ipart];
00145       // and merge them (with a fake dij)
00146       double dij = 0.0;
00147 
00148       // create the new jet by hand so that we can adjust its user index
00149       PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
00150 
00151       // set the user index to be the pass in which the jet was discovered
00152       newjet.set_user_index(jet.pass);
00153         
00154       clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k);
00155     }
00156     // we have merged all the jet's particles into a single object, so now
00157     // "declare" it to be a beam (inclusive) jet.
00158     // [NB: put a sensible looking d_iB just to be nice...]
00159     double d_iB = clust_seq.jets()[jet_k].perp2();
00160     clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
00161   }
00162 
00163   // now copy the list of protocones into an "extras" objects
00164   SISConeExtras * extras = new SISConeExtras;
00165   for (unsigned ipass = 0; ipass < siscone->protocones_list.size(); ipass++) {
00166     for (unsigned ipc = 0; ipc < siscone->protocones_list[ipass].size(); ipc++) {
00167       double rap = siscone->protocones_list[ipass][ipc].eta;
00168       double phi = siscone->protocones_list[ipass][ipc].phi;
00169       PseudoJet protocone(cos(phi),sin(phi),sinh(rap),cosh(rap));
00170       //PseudoJet protocone(siscone->protocones_list[ipass][ipc]);
00171       protocone.set_user_index(ipass);
00172       extras->_protocones.push_back(protocone);
00173     }
00174   }
00175   extras->_most_ambiguous_split = siscone->most_ambiguous_split;
00176 
00177   // tell it what the jet definition was
00178   extras->_jet_def_plugin = this;
00179 
00180   // give the extras object to the cluster sequence.
00181   clust_seq.plugin_associate_extras(auto_ptr<ClusterSequence::Extras>(extras));
00182 }
00183 
00184 
00186 string SISConeExtras::description() const {
00187   ostringstream ostr;
00188   ostr << "This SISCone clustering found " << protocones().size() 
00189        << " stable protocones";
00190   return ostr.str();
00191 }
00192 
00193 // OBSOLETE CODE  
00195 //if (siscone.n_warnings > 0) {
00196 //  // print the event (very dirty...)
00197 //  cout.precision(14);
00198 //  cout << "Culprit event is:" << endl;
00199 //  for (int i = 0; i < n; i++) {
00200 //    const PseudoJet & p = clust_seq.jets()[i]; // shorthand
00201 //    cout << p.px() << " " << p.py() << " " << p.pz() << " " << p.E() << endl;
00202 //  }
00203 //}
00204 
00205 
00206 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh

Generated on Tue Dec 18 17:05:02 2007 for fastjet by  doxygen 1.5.2