source: svn/trunk/Delphes.cpp@ 36

/*
  ---- Delphes ----
  A Fast Simulator for general purpose LHC detector
  S. Ovyn ~~~~ severine.ovyn@uclouvain.be

  Center for Particle Physics and Phenomenology (CP3)
  Universite Catholique de Louvain (UCL)
  Louvain-la-Neuve, Belgium
*/

/// \file Delphes.cpp
/// \brief executable for the Delphes

#include "TChain.h"
#include "TApplication.h"

#include "Utilities/ExRootAnalysis/interface/ExRootTreeReader.h"
#include "Utilities/ExRootAnalysis/interface/ExRootTreeWriter.h"
#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"

#include "H_BeamParticle.h"
#include "H_BeamLine.h"
#include "H_RomanPot.h"

#include "interface/DataConverter.h"
#include "interface/HEPEVTConverter.h"
#include "interface/LHEFConverter.h"
#include "interface/STDHEPConverter.h"

#include "interface/SmearUtil.h"
#include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
#include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh"

// get info on how fastjet was configured
#include "Utilities/Fastjet/include/fastjet/config.h"

// make sure we have what is needed
#ifdef ENABLE_PLUGIN_SISCONE
#  include "Utilities/Fastjet/plugins/SISCone/SISConePlugin.hh"
#endif
#ifdef ENABLE_PLUGIN_CDFCONES
#  include "Utilities/Fastjet/plugins/CDFCones/CDFMidPointPlugin.hh"
#  include "Utilities/Fastjet/plugins/CDFCones/CDFJetCluPlugin.hh"
#endif

#include<vector>
#include<iostream>



using namespace std;

//------------------------------------------------------------------------------
void todo(string filename) {
  ifstream infile(filename.c_str());
  cout << "** TODO list ..." << endl;
  while(infile.good()) {
        string temp;
        getline(infile,temp);
        cout << "*" << temp << endl;
  }
  cout << "** done...\n"; 
}

//------------------------------------------------------------------------------

int main(int argc, char *argv[])
{
  int appargc = 2;
  char *appName = "JetsSim";
  char *appargv[] = {appName, "-b"};
  TApplication app(appName, &appargc, appargv);
  
  if(argc != 4 && argc != 3) {
      cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
      cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
      cout << " output_file - output file." << endl;
      cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
      exit(1);
  }

  srand (time (NULL));         /* Initialisation du générateur */
  
  //read the input TROOT file
  string inputFileList(argv[1]), outputfilename(argv[2]);
  if(outputfilename.find(".root") > outputfilename.length() ) {
        cout << "output_file should be a .root file!\n"; 
        exit(1);
  }

  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
  outputFile->Close();

  string line;
  ifstream infile(inputFileList.c_str());
  infile >> line; // the first line determines the type of input files
 
  DataConverter *converter=0;
   
  if(strstr(line.c_str(),".hep"))
    {                           
      cout<<"*************************************************************************"<<endl;
      cout<<"************         StdHEP file format detected           **************"<<endl;
      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
      cout<<"*************************************************************************"<<endl;
      converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
    }
  else if(strstr(line.c_str(),".lhe"))
    {
      cout<<"*************************************************************************"<<endl;
      cout<<"************          LHEF file format detected            **************"<<endl;
      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
      cout<<"*************************************************************************"<<endl;
      converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
    }
  else if(strstr(line.c_str(),".root"))
    {
      cout<<"*************************************************************************"<<endl;
      cout<<"************         h2root file format detected           **************"<<endl;
      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
      cout<<"*************************************************************************"<<endl;
      converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
    }
  else { cout << "***  " << line.c_str() << "\n***  file format not identified\n***  Exiting\n"; return -1;};
  
  TChain chain("GEN");
  chain.Add(outputfilename.c_str());
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  const TClonesArray *branchGen = treeReader->UseBranch("Particle");
  TIter itGen((TCollection*)branchGen);
  
  //write the output root file
  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
  ExRootTreeBranch *branchJet = treeWriter->NewBranch("Jet", TRootJet::Class());
  ExRootTreeBranch *branchTauJet = treeWriter->NewBranch("TauJet", TRootTauJet::Class());
  ExRootTreeBranch *branchElectron = treeWriter->NewBranch("Electron", TRootElectron::Class());
  ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class());
  ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class());
  ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::Class());
  ExRootTreeBranch *branchETmis = treeWriter->NewBranch("ETmis", TRootETmis::Class());
  ExRootTreeBranch *branchCalo = treeWriter->NewBranch("CaloTower", TRootCalo::Class());
  ExRootTreeBranch *branchZDC = treeWriter->NewBranch("ZDChits", TRootZdcHits::Class());
  ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
  ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
  
  
  TRootGenParticle *particle; 
  TRootETmis *elementEtmis;
  TRootElectron *elementElec;
  TRootMuon *elementMu;
  TRootPhoton *elementPhoton;
  TRootJet *elementJet;
  TRootTauJet *elementTauJet;
  TRootTracks *elementTracks;
  TRootCalo *elementCalo;
  TRootZdcHits *elementZdc;
  TRootRomanPotHits *elementRP220, *elementFP420;
  
  //read the datacard input file
  string DetDatacard("");
  if(argc==4)  DetDatacard =argv[3];
  RESOLution *DET = new RESOLution();   
  DET->ReadDataCard(DetDatacard);
  
  TLorentzVector genMomentum(0,0,0,0);
  LorentzVector jetMomentum;
  vector<TLorentzVector> TrackCentral;  
  vector<PhysicsTower> towers;
  
  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
  vector<fastjet::PseudoJet> inclusive_jets;
  vector<fastjet::PseudoJet> sorted_jets;
  
  vector<TLorentzVector> electron;
  vector<int> elecPID;
  vector<TLorentzVector> muon;
  vector<int> muonPID;
  TSimpleArray<TRootGenParticle> NFCentralQ;
  
  //Initialisation of Hector
  extern bool relative_energy;
  relative_energy = true; // should always be true
  extern int kickers_on;
  kickers_on = 1;         // should always be 1
  
  // user should provide : (1) optics file for each beamline, and IPname, 
  // and offset data (s,x) for optical elements
  H_BeamLine* beamline1 = new H_BeamLine(1,500.);
  beamline1->fill("data/LHCB1IR5_v6.500.tfs",1,"IP5");
  beamline1->offsetElements(120,-0.097);
  H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",220,2000); // RP 220m, 2mm, beam 1
  H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",420,4000); // RP 420m, 4mm, beam 1
  beamline1->add(rp220_1);
  beamline1->add(rp420_1);
  
  H_BeamLine* beamline2 = new H_BeamLine(1,500.);
  beamline2->fill("data/LHCB1IR5_v6.500.tfs",-1,"IP5");
  beamline2->offsetElements(120,+0.097);
  H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",220,2000);// RP 220m, 2mm, beam 2
  H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",420,4000);// RP 420m, 4mm, beam 2
  beamline2->add(rp220_2);
  beamline2->add(rp420_2);
  
  // we will have four jet definitions, and the first three will be
  // plugins
  fastjet::JetDefinition jet_def;
  fastjet::JetDefinition::Plugin * plugins;
  
  switch(DET->JETALGO) {
  default:
  case 1: {
    
    // set up a CDF midpoint jet definition
#ifdef ENABLE_PLUGIN_CDFCONES
    plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    jet_def = fastjet::JetDefinition(plugins);
#else
    plugins = NULL;
#endif
  }
    break;
    
  case 2: {
    
    // set up a CDF midpoint jet definition
#ifdef ENABLE_PLUGIN_CDFCONES
    plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
    jet_def = fastjet::JetDefinition(plugins);
#else
    plugins = NULL;
#endif
  }
    break;
  case 3: {
    // set up a siscone jet definition
#ifdef ENABLE_PLUGIN_SISCONE
    int npass = 0;               // do infinite number of passes
    double protojet_ptmin = 0.0; // use all protojets
    plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
    jet_def = fastjet::JetDefinition(plugins);
#else
    plugins = NULL;
#endif
  }
    break;
    
  case 4: {
    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
    //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
    //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
  }
    break;
  }
  
  // Loop over all events
  Long64_t entry, allEntries = treeReader->GetEntries();
  cout << "** Chain contains " << allEntries << " events" << endl;
  for(entry = 0; entry < allEntries; ++entry)
    {
      TLorentzVector PTmis(0,0,0,0);
      treeReader->ReadEntry(entry);
      treeWriter->Clear();
      
      if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
      
      electron.clear();
      muon.clear();
      elecPID.clear();
      muonPID.clear();
      NFCentralQ.Clear();
      
      itGen.Reset();
      TrackCentral.clear();
      towers.clear();
      input_particles.clear();
      inclusive_jets.clear();
      sorted_jets.clear();
      
      // Loop over all particles in event
      while( (particle = (TRootGenParticle*) itGen.Next()) )
        {
          genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
          
          int pid = abs(particle->PID);
          float eta = fabs(particle->Eta);
          
          //subarray of the quarks (i.e. not final particles, i.e status not equal to 1) 
          // in the tracker (i.e. for those we know the tracks)
          // with enough p_T
          //// This subarray is needed for the B-jet algorithm
          // optimization for speed : put first PID condition, then ETA condition, then either pt or status
          if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ?
              eta < DET->MAX_TRACKER &&
              particle->Status != 1 &&
              particle->PT > DET->PT_QUARKS_MIN ) {
            NFCentralQ.Add(particle);
          }
          
          
          // keeps only final particles, visible by the central detector, including the fiducial volume
          // the ordering of conditions have been optimised for speed : put first the STATUS condition
          if( (particle->Status == 1)    &&
              ( 
               (pid == pMU  && eta < DET->MAX_MU) ||
               (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) 
                )
              ) { 
            switch(pid) {
              
            case pE: // all electrons with eta < DET->MAX_CALO_FWD
              DET->SmearElectron(genMomentum);
              electron.push_back(genMomentum);
              elecPID.push_back(particle->PID);
              break; // case pE
            case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
              DET->SmearElectron(genMomentum);
              if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
                elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
                elementPhoton->Set(genMomentum);
              }
              break; // case pGAMMA
            case pMU: // all muons with eta < DET->MAX_MU
              DET->SmearMu(genMomentum);
              muonPID.push_back(particle->PID);
              muon.push_back(genMomentum);
              break; // case pMU
            case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
            case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
              DET->SmearHadron(genMomentum, 0.7);
              break; // case hadron
            default:   // all other final particles with eta < DET->MAX_CALO_FWD
              DET->SmearHadron(genMomentum, 1.0);
              break;
            } // switch (pid)
            
            // all final particles but muons and neutrinos      
            // for calorimetric towers and mission PT
            if(genMomentum.E()!=0) {
              PTmis = PTmis + genMomentum;//ptmis
              if(pid !=pMU) {
                
                towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
                // create a fastjet::PseudoJet with these components and put it onto
                // back of the input_particles vector
                input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
                elementCalo = (TRootCalo*) branchCalo->NewEntry();
                elementCalo->Set(genMomentum);
              }
            }
            
            
            // all final charged particles
            if(
               ((rand()%100) < DET->TRACKING_EFF)  &&
               (genMomentum.E()!=0) &&
               (fabs(particle->Eta) < DET->MAX_TRACKER) && 
               (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
               (pid != pGAMMA)    && 
               (pid != pPI0)      &&
               (pid != pK0L)      &&
               (pid != pN)        &&
               (pid != pSIGMA0)   &&
               (pid != pDELTA0)   &&
               (pid != pK0S)   // not charged particles : invisible by tracker
               )
              {
                elementTracks = (TRootTracks*) branchTracks->NewEntry();
                elementTracks->Set(genMomentum);
                TrackCentral.push_back(genMomentum);
              }
          } // switch
          
          // Forward particles in CASTOR ?
          /*      if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) 
                  && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
                  
                  
                  } // CASTOR
          */
          // Zero degree calorimeter, for forward neutrons and photons
          if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > DET->MIN_ZDC ) {
            // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
            // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
            // !!!!!!!!! pas de smearing ! on garde trop d'info !
            elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
            elementZdc->Set(genMomentum);
            
            // time of flight t is t = T + d/[ cos(theta) v ]
            //double tx = acos(particle->Px/particle->Pz);
            //double ty = acos(particle->Py/particle->Pz);
            //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
            //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
            double flight_distance = (DET->ZDC_S - particle->Z*(1E-6)); 
            // assumes also that the emission angle is so small that 1/(cos theta) = 1
            elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles 
            elementZdc->side = sign(particle->Eta);
            
          }
          
          // if forward proton
          if( (pid == pP)  && (particle->Status == 1) &&  (fabs(particle->Eta) > DET->MAX_CALO_FWD) ) 
            {
              // !!!!!!!! put here particle->CHARGE and particle->MASS
              H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
              p1.smearAng();
              p1.smearPos();
              p1.setPosition(p1.getX()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,p1.getTY(),0);
              p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
              
              H_BeamLine *beamline;
              if(particle->Eta >0) beamline = beamline1;
              else beamline = beamline2;
              
              p1.computePath(beamline,1);
              
              if(p1.stopped(beamline)) {
                if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
                  p1.propagate(DET->RP220_S);
                  elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
                  elementRP220->X  = (1E-6)*p1.getX();  // [m]
                  elementRP220->Y  = (1E-6)*p1.getY();  // [m]
                  elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
                  elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
                  elementRP220->S = p1.getS();          // [m]
                  elementRP220->T = -1;                 // not yet implemented
                  elementRP220->E = p1.getE();  // not yet implemented
                  elementRP220->q2 = -1;                // not yet implemented
                  elementRP220->side = sign(particle->Eta); 
                  
                } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
                  p1.propagate(DET->FP420_S);
                  elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
                  elementFP420->X  = (1E-6)*p1.getX();  // [m]
                  elementFP420->Y  = (1E-6)*p1.getY();  // [m]
                  elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
                  elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
                  elementFP420->S = p1.getS();        // [m]
                  elementFP420->T = -1;                       // not yet implemented
                  elementFP420->E = p1.getE();          // not yet implemented
                  elementFP420->q2 = -1;                // not yet implemented
                  elementFP420->side = sign(particle->Eta);
                }
              }
              
              //                if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
              //                        cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
            } // if forward proton
          
        } // while
      
      for(unsigned int i=0; i < electron.size(); i++) {
        if(electron[i].E()!=0 && fabs(electron[i].Eta()) < DET->MAX_TRACKER && electron[i].Pt() > DET->ELEC_pt)
          {
            elementElec = (TRootElectron*) branchElectron->NewEntry();
            elementElec->Set(electron[i]);
            elementElec->Charge = sign(elecPID[i]);
            elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);
          }
      }
      for(unsigned int i=0; i < muon.size(); i++) {
        if(muon[i].E()!=0 && fabs(muon[i].Eta()) < DET->MAX_MU && muon[i].Pt() > DET->MUON_pt) 
          {
            elementMu = (TRootMuon*) branchMuon->NewEntry();
            elementMu->Charge = sign(muonPID[i]);
            elementMu->Set(muon[i]);
            elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,2.0);
          }
      }
      
      
      // computes the Missing Transverse Momentum
      elementEtmis = (TRootETmis*) branchETmis->NewEntry();
      elementEtmis->ET = (-PTmis).Pt();
      elementEtmis->Phi = (-PTmis).Phi();
      elementEtmis->Px = (-PTmis).Px();
      elementEtmis->Py = (-PTmis).Py();
      
      //*****************************
      
      // run the jet clustering with the above jet definition
      if(input_particles.size()!=0)
        {
          fastjet::ClusterSequence clust_seq(input_particles, jet_def);
          
          
          // extract the inclusive jets with pt > 5 GeV
          double ptmin = 5.0;
          inclusive_jets = clust_seq.inclusive_jets(ptmin);
          
          // sort jets into increasing pt
          sorted_jets = sorted_by_pt(inclusive_jets);
        }
      for (unsigned int i = 0; i < sorted_jets.size(); i++) {
        elementJet = (TRootJet*) branchJet->NewEntry();
        TLorentzVector JET;
        JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
        // Tau jet identification : 1! track and electromagnetic collimation
        if(fabs(JET.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) {
        double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
        if(
           ( Energie_tau_central/JET.E() > DET->TAU_EM_COLLIMATION ) &&
           ( DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 ) &&
           ( JET.Pt() > DET->TAUJET_pt)
          ) {
               elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
               elementTauJet->Set(JET);
            } // if tau jet
           } // if JET.eta < tracker - tau_cone : Tau jet identification
  
        if(JET.Pt() > DET->JET_pt)
          {
            elementJet->Set(JET);
            // b-jets
            bool btag=false;
            if((fabs(JET.Eta()) < DET->MAX_TRACKER && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
            elementJet->Btag = btag;
          }
      } // for itJet : loop on all jets
      
      treeWriter->Fill();
      // Add here the trigger
      // Should test all the trigger table on the event, based on reconstructed objects
    } // Loop over all events
  treeWriter->Write();
  
  cout << "** Exiting..." << endl;
  
  delete treeWriter;
  delete treeReader;
  delete DET;
  if(converter) delete converter;

  todo("TODO");
}