stableTest2.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"

#include "TFile.h"
#include "TTree.h"


//-------- 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] dataDir" << 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;
    } 

  char * dataDir = 0;
 
  for ( int index = optind; index < argc; index++ ) {
    if ( index == optind ) dataDir = argv[index];
  }

  if ( !dataDir ) {
    std::cerr << "Must provide data directory as non-option argument" << std::endl;
    return 1;
  }  

  // -------------- root tree branch variables
  TFile f("output.root","RECREATE");
  TTree *tree = new TTree("tree","tree");
  
  std::vector<double> *pt=0;
  std::vector<int> *pid=0;

  double bpt,bphi,beta,bbarpt,bbarphi,bbareta;
  int bstatus,bbarstatus;
  double bdR;
  
  tree->Branch("pt",&pt);
  tree->Branch("pid",&pid);
  tree->Branch("bpt",&bpt,"bpt/D");
  tree->Branch("bphi",&bphi,"bphi/D");
  tree->Branch("beta",&beta,"beta/D");
  tree->Branch("bstatus",&bstatus,"bstatus/I");
  tree->Branch("bbarpt",&bbarpt,"bbarpt/D");
  tree->Branch("bbarphi",&bbarphi,"bbarphi/D");
  tree->Branch("bbareta",&bbareta,"bbareta/D");
  tree->Branch("bbarstatus",&bbarstatus,"bbarstatus/I");
  tree->Branch("bdR",&bdR,"bdR/D");


  //---------------- 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++ ) {

    pt->clear();
    pid->clear();

    stab::Event event;
    stab::StabilityEventInfo evInfo;
    
    // restore event evNum from archive
    ref.restore(evNum, &ev); 
    SimpleEvent & evs = ev.sim;
    const unsigned nSimulatedJets = evs.size();

    PythiaEventInfo & pyInfo = ev.pyev;
    const int pyN = pyInfo.N();
    const int bnum = pyInfo.findFirst(5);
    const int bbarnum = pyInfo.findFirst(-5);

    bpt = pyInfo.pt(bnum);
    bphi = atan2(pyInfo.py(bnum),pyInfo.px(bnum));
    beta = asinh(pyInfo.pz(bnum)/bpt);
    bstatus = pyInfo.statusCode(bnum);

    bbarpt = pyInfo.pt(bbarnum);
    bbarphi = atan2(pyInfo.py(bbarnum),pyInfo.px(bbarnum));
    bbareta = asinh(pyInfo.pz(bbarnum)/bbarpt);
    bbarstatus = pyInfo.statusCode(bbarnum);

    const double dEta = beta - bbareta;
    double dPhi = bphi - bbarphi;
    while ( dPhi > M_PI) dPhi -= 2*M_PI;
    while ( dPhi <= -M_PI ) dPhi += 2*M_PI; 

    bdR = sqrt(dEta*dEta+dPhi*dPhi);
 
    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 );
        pt->push_back(pj.pt());
    pid->push_back(part.code());
      }
    }

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

    tree->Fill();
  }

  /*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;
  }*/

  tree->Write();
  f.Close();
 
  //delete tree; 
    
return 0;
}