stableTest.cc

//a program to analyze jet stability


#include <iostream>
#include <fstream>
#include <cstdlib>
#include <unistd.h>
#include <vector>
#include <algorithm>

#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"

#include "JetStabilityDefs.h"
#include "JetStability.h"
#include "DiffusionDistFastJet.h"
#include "SimpleSmear.h"
#include "SimpleEvent.hh"
#include "SimulatedEvent.hh"
#include "QCDandCmsLikeSmear.h"

#include "PythiaJetGun.hh"

#include "geners/BinaryFileArchive.hh"
#include "geners/Reference.hh"


//-------- teypdefs and stuff -------------------
using namespace gs;

bool oddSort ( const std::pair<unsigned,double> & a, const std::pair<unsigned,double> &b ) {
  return a.second > b.second;
}


void printUsage(int argc, char **argv ) {
  if ( argv[0] )
    std::cout << argv[0] << " -a archive [-A nEvents] [-C nCycles] [-o outputFile] [-q fuzzyFile] [-m scanMin] [-M scanMax] [-n nScan] [-f fixedVar] [-E (set exclusive)] [-d distCut] [-N nJets] [-r R]" << std::endl;
}


int main (int argc, char **argv ) {


  // options loaded from command line options
  char *archiveName = 0;
  unsigned nEvents = 0;
  char *outFile = 0;
  char *fuzzyFile = 0;
  bool exclusive = false;
  double distCut=1.;
  unsigned N = 2U;
  double R = 0.5;
  unsigned nCycles=100;

  double scanMin = 0.1;
  double scanMax = 1.;
  unsigned scanN = 10;

  double fixed = 1.;
  
  //double qcdCutOff = exp(3.);
  //double qcdCutOff = 0.01;
  double qcdCutOff = 0.705796; // 75% quantile of dist for FluctuateQCD2_geom
  bool doFlatPhaseSpace = true;
 
  int c;

  while (( c = getopt(argc,argv,"a:A:o:r:C:N:m:M:n:f:Ed:q:")) != -1 )
    switch (c) {
      case 'a':
    archiveName = optarg;
    break;
      case 'A':
    nEvents = atoi(optarg);
    break;
      case 'N':
    N = atoi(optarg);
    break;
      case 'r':
    R = atof(optarg);
    break;
      case 'o':
    outFile = optarg;
    break;
      case 'C':
    nCycles = atoi(optarg);
    break;
      case 'm':
    scanMin = atof(optarg);
    break;
      case 'M':
    scanMax = atof(optarg);
    break;
      case 'n':
    scanN = atoi(optarg);
    break;
      case 'f':
    fixed = atof(optarg);
    break;
      case 'E':
    exclusive = true;
    break;
      case 'd':
    distCut = atof(optarg);
    break;
      case 'q':
    fuzzyFile = optarg;
    break;
      default:
    printUsage(argc,argv);
    return 0;
    } 

  //---------------- the program ----------------------------

  // stability output stream
  std::ofstream out;
  if ( outFile ) 
    out.open(outFile,std::ios_base::out);

  // fuzzyness output stream
  std::ofstream fout;
  if ( fuzzyFile ) 
    fout.open(fuzzyFile,std::ios_base::out);
   

  std::cout << "Opening archive.......";
  // Open the archive
  BinaryFileArchive ar(archiveName, "r");
  if (!ar.isOpen())
  {
      std::cerr << ar.error() << std::endl;
      return 1;
  }
  std::cout << "done" << std::endl;

  // Find all events in the archive with the given category and process them
  std::cout << "Finding events.......";
  Reference<SimulatedEvent> ref(ar, Regex(".*"), "dijet");
  if ( !nEvents ) nEvents = ref.size();
  std::cout << "done  found: " << ref.size() << " events" << std::endl;

  // fastjet jet definition
  stab::DiffusionDistFastJet fjDiffDist(R);
  stab::JetStability jstab((fastjet::JetDefinition*)&fjDiffDist);
  jstab.setScan(scanN,scanMin,scanMax);
  jstab.setNCycles(nCycles);

  // simple smearing agent
  stab::SimpleSmear smear(0.087,0.15);
  //stab::QCDandCmsLikeSmear smear;
  

  SimulatedEvent ev;
  // loop over so many events
  for ( unsigned evNum=0; evNum != nEvents; evNum++ ) {

    stab::Event event;
    stab::StabilityEventInfo evInfo;
 
    // restore event evNum from archive
    ref.restore(evNum, &ev); 
    SimpleEvent & evs = ev.sim;
    const unsigned nSimulatedJets = evs.size();
    for ( unsigned j=0; j != nSimulatedJets; j++ ) {
      const std::vector<PythiaJetGun::Particle> & constituents = evs[j].jetParticles;
      const unsigned nConstituents = constituents.size();
      for ( unsigned c=0; c != nConstituents; c++ ) {
    const PythiaJetGun::Particle & part = constituents[c];
        const rk::P4 & p4 = part.p4(); 
        const geom3::Vector3 & mom = p4.momentum();
    stab::PseudoJet pj(mom.x(),mom.y(),mom.z(),p4.e());
    //stab::StabilityUserInfo *pjInfo = new stab::StabilityUserInfo(part.code(),part.charge());
        //pj.set_user_info( pjInfo );
    evInfo.push_back( stab::StabilityUserInfo(part.code(),part.charge()) ); 
    event.push_back( pj );
      }
    }

    jstab.reset(event,evInfo); 
    jstab.run(smear); 

  }

  const std::vector<std::vector<double> > & stability = jstab.pairwiseRandStability();
  const std::vector<std::vector<double> > & fuzziness = jstab.fuzzinessNJet();

  const double sMax = jstab.scanMax();
  const double sMin = jstab.scanMin();
  const unsigned sN = jstab.nScan();
  const double dScan = (sMax-sMin)/sN;

  const unsigned nStab = stability.size();
  assert( nStab == sN );
  assert( stability.size() == fuzziness.size() );
  for ( unsigned s=0; s != nStab; s++ ) {
    const double dCut = sMin + dScan*s;
    const std::vector<double> & stabPoints = stability[s];
    const std::vector<double> & fuzzPoints = fuzziness[s];
    const unsigned nPoints = stabPoints.size();
    long double avg = 0.L;
    long double favg = 0.L;
    for ( unsigned p=0; p != nPoints; p++ ) {
      avg += stabPoints[p];
      favg += fuzzPoints[p];  
    }
    avg /= nPoints;
    favg /= nPoints;
    out << dCut << " " << avg << std::endl;
    fout << dCut << " " << favg << std::endl;
  } 
    
return 0;
}