Fork me on GitHub

Changeset 19 in svn for trunk


Ignore:
Timestamp:
Nov 6, 2008, 5:48:51 PM (16 years ago)
Author:
severine ovyn
Message:

plugins at the begining

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/Delphes.cpp

    r16 r19  
    197197  beamline2->add(rp420_2);
    198198
    199  
     199  // we will have four jet definitions, and the first three will be
     200  // plugins
     201  fastjet::JetDefinition jet_def;
     202  fastjet::JetDefinition::Plugin * plugins;
     203
     204  switch(DET->JETALGO) {
     205  default:
     206  case 1: {
     207
     208   // set up a CDF midpoint jet definition
     209   #ifdef ENABLE_PLUGIN_CDFCONES
     210   plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
     211   jet_def = fastjet::JetDefinition(plugins);
     212   #else
     213   plugins = NULL;
     214   #endif
     215   }
     216   break;
     217
     218   case 2: {
     219
     220   // set up a CDF midpoint jet definition
     221   #ifdef ENABLE_PLUGIN_CDFCONES
     222   plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
     223   jet_def = fastjet::JetDefinition(plugins);
     224   #else
     225   plugins = NULL;
     226   #endif
     227   }
     228   break;
     229   case 3: {
     230   // set up a siscone jet definition
     231   #ifdef ENABLE_PLUGIN_SISCONE
     232   int npass = 0;               // do infinite number of passes
     233   double protojet_ptmin = 0.0; // use all protojets
     234   plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
     235   jet_def = fastjet::JetDefinition(plugins);
     236   #else
     237   plugins = NULL;
     238   #endif
     239   }
     240   break;
     241
     242  case 4: {
     243  jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
     244  //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
     245  //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
     246  }
     247  break;
     248  }
    200249 
    201   // Loop over all events
     250   // Loop over all events
    202251  Long64_t entry, allEntries = treeReader->GetEntries();
    203252  cout << "** Chain contains " << allEntries << " events" << endl;
     
    409458      //*****************************
    410459
    411       // we will have four jet definitions, and the first three will be
    412       // plugins
    413       fastjet::JetDefinition jet_def;
    414       fastjet::JetDefinition::Plugin * plugins;
    415      
    416       switch(DET->JETALGO) {
    417       default:
    418       case 1: {
    419        
    420         // set up a CDF midpoint jet definition
    421         #ifdef ENABLE_PLUGIN_CDFCONES
    422         plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    423         jet_def = fastjet::JetDefinition(plugins);
    424         #else
    425         plugins = NULL;
    426         #endif
    427        
    428       }
    429         break;
    430        
    431       case 2: {
    432        
    433         // set up a CDF midpoint jet definition
    434         #ifdef ENABLE_PLUGIN_CDFCONES
    435         plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
    436         jet_def = fastjet::JetDefinition(plugins);
    437         #else
    438         plugins = NULL;
    439         #endif
    440       }
    441         break;
    442        
    443       case 3: {
    444         // set up a siscone jet definition
    445         #ifdef ENABLE_PLUGIN_SISCONE
    446         int npass = 0;               // do infinite number of passes
    447         double protojet_ptmin = 0.0; // use all protojets
    448         plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
    449         jet_def = fastjet::JetDefinition(plugins);
    450         #else
    451         plugins = NULL;
    452         #endif
    453       }
    454         break;
    455        
    456       case 4: {
    457         jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
    458         //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
    459         //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
    460       }
    461         break;
    462       }
    463460      // run the jet clustering with the above jet definition
    464461      if(input_particles.size()!=0)
     
    478475        TLorentzVector JET;
    479476        JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
     477        cout<<"Jet.Pt() "<<JET.Pt()<<endl;
    480478        elementJet->Set(JET);
    481479        // b-jets
  • trunk/Makefile

    r15 r19  
    5656        interface/LHEFConverter.h \
    5757        interface/STDHEPConverter.h \
    58         interface/SmearUtil.h
     58        interface/SmearUtil.h \
     59        Utilities/Fastjet/include/fastjet/PseudoJet.hh \
     60        Utilities/Fastjet/include/fastjet/ClusterSequence.hh \
     61        Utilities/Fastjet/include/fastjet/config.h \
     62        interface/TreeClasses.h
    5963EXECUTABLE =  \
    6064        Delphes$(ExeSuf) \
     
    6872        Utilities/ExRootAnalysis/src/BlockClassesLinkDef.h \
    6973        Utilities/ExRootAnalysis/interface/BlockClasses.h
     74tmp/src/TreeClassesDict.$(SrcSuf): \
     75        src/TreeClassesLinkDef.h \
     76        interface/TreeClasses.h
    7077DICT =  \
    71         tmp/Utilities/ExRootAnalysis/src/BlockClassesDict.$(SrcSuf)
     78        tmp/Utilities/ExRootAnalysis/src/BlockClassesDict.$(SrcSuf) \
     79        tmp/src/TreeClassesDict.$(SrcSuf)
    7280
    7381DICT_OBJ =  \
    74         tmp/Utilities/ExRootAnalysis/src/BlockClassesDict.$(ObjSuf)
     82        tmp/Utilities/ExRootAnalysis/src/BlockClassesDict.$(ObjSuf) \
     83        tmp/src/TreeClassesDict.$(ObjSuf)
    7584
    7685tmp/src/HEPEVTConverter.$(ObjSuf): \
     
    101110        src/SmearUtil.$(SrcSuf) \
    102111        interface/SmearUtil.h
     112tmp/src/TreeClasses.$(ObjSuf): \
     113        src/TreeClasses.$(SrcSuf) \
     114        interface/TreeClasses.h
    103115tmp/Utilities/ExRootAnalysis/src/BlockClasses.$(ObjSuf): \
    104116        Utilities/ExRootAnalysis/src/BlockClasses.$(SrcSuf) \
     
    281293        tmp/src/STDHEPConverter.$(ObjSuf) \
    282294        tmp/src/SmearUtil.$(ObjSuf) \
     295        tmp/src/TreeClasses.$(ObjSuf) \
    283296        tmp/Utilities/ExRootAnalysis/src/BlockClasses.$(ObjSuf) \
    284297        tmp/Utilities/ExRootAnalysis/src/ExRootProgressBar.$(ObjSuf) \
     
    468481        @touch $@
    469482
     483interface/TreeClasses.h: \
     484        Utilities/ExRootAnalysis/interface/BlockCompare.h \
     485        Utilities/ExRootAnalysis/interface/BlockClasses.h
     486        @touch $@
     487
    470488interface/LHEFConverter.h: \
    471489        Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h \
     
    497515        @touch $@
    498516
     517Utilities/ExRootAnalysis/interface/BlockClasses.h: \
     518        Utilities/ExRootAnalysis/interface/BlockCompare.h
     519        @touch $@
     520
    499521interface/STDHEPConverter.h: \
    500522        Utilities/ExRootAnalysis/interface/BlockClasses.h \
     
    502524        Utilities/ExRootAnalysis/interface/LHEF.h \
    503525        interface/DataConverter.h
    504         @touch $@
    505 
    506 Utilities/ExRootAnalysis/interface/BlockClasses.h: \
    507         Utilities/ExRootAnalysis/interface/BlockCompare.h
    508526        @touch $@
    509527
  • trunk/Resolutions.cpp

    r9 r19  
    1919#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
    2020
    21 #include "Utilities/CDFCones/interface/JetCluAlgorithm.h"
    22 #include "Utilities/CDFCones/interface/MidPointAlgorithm.h"
    23 #include "Utilities/CDFCones/interface/PhysicsTower.h"
    24 #include "Utilities/CDFCones/interface/Cluster.h"
     21#include "H_BeamParticle.h"
     22#include "H_BeamLine.h"
     23#include "H_RomanPot.h"
    2524
    2625#include "interface/DataConverter.h"
     
    2827#include "interface/LHEFConverter.h"
    2928#include "interface/STDHEPConverter.h"
     29
    3030#include "interface/SmearUtil.h"
     31#include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
     32#include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh"
     33
     34// get info on how fastjet was configured
     35#include "Utilities/Fastjet/include/fastjet/config.h"
     36
     37// make sure we have what is needed
     38#ifdef ENABLE_PLUGIN_SISCONE
     39#  include "Utilities/Fastjet/plugins/SISCone/SISConePlugin.hh"
     40#endif
     41#ifdef ENABLE_PLUGIN_CDFCONES
     42#  include "Utilities/Fastjet/plugins/CDFCones/CDFMidPointPlugin.hh"
     43#  include "Utilities/Fastjet/plugins/CDFCones/CDFJetCluPlugin.hh"
     44#endif
     45
     46#include<vector>
     47#include<iostream>
     48
    3149#include "interface/TreeClasses.h"
    32 
    3350using namespace std;
    3451
     
    3754
    3855
    39 void PairingJet(TLorentzVector &JETSm, TLorentzVector JET, vector<Cluster> jetsS)
     56void PairingJet(TLorentzVector &JETSm, TLorentzVector JET,  vector<fastjet::PseudoJet> sorted_jetsS)
    4057{
    4158  JETSm.SetPxPyPzE(0,0,0,0);
    42   vector<Cluster>::iterator itJetS;
    4359  float deltaRtest=5000;
    44   LorentzVector jetMomentumS;
    45   for(itJetS = jetsS.begin(); itJetS != jetsS.end(); ++itJetS)
    46      {
    47        jetMomentumS = itJetS->fourVector;
    48        TLorentzVector Att;
    49        Att.SetPxPyPzE(jetMomentumS.px,jetMomentumS.py,jetMomentumS.pz,jetMomentumS.E);
    50        if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
    51             {
    52               deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
    53               if(deltaRtest < 0.25)
    54                 {
    55                   JETSm = Att;
    56                 }
    57             }
    58      }
     60  for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
     61      TLorentzVector Att;
     62      Att.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
     63      if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
     64        {
     65          deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
     66          if(deltaRtest < 0.25)
     67             {
     68               JETSm = Att;
     69             }
     70        }
     71    }
    5972}
    6073
     
    148161  DET->ReadDataCard(DetDatacard);
    149162 
    150  
    151163  TLorentzVector genMomentum(0,0,0,0);
    152164  LorentzVector jetMomentum;
    153   vector<TLorentzVector> TrackCentral; 
    154   vector<PhysicsTower> towersS;
    155   vector<Cluster> jetsS;
    156   vector<PhysicsTower> towers;
    157   vector<Cluster> jets;
     165  vector<TLorentzVector> TrackCentral;
     166
     167  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
     168  vector<fastjet::PseudoJet> inclusive_jets;
     169  vector<fastjet::PseudoJet> sorted_jets;
     170
     171  vector<fastjet::PseudoJet> input_particlesS;//for FastJet algorithm
     172  vector<fastjet::PseudoJet> inclusive_jetsS;
     173  vector<fastjet::PseudoJet> sorted_jetsS;
     174
     175  fastjet::JetDefinition jet_def;
     176  fastjet::JetDefinition jet_defS;
     177  fastjet::JetDefinition::Plugin * plugins;
     178  fastjet::JetDefinition::Plugin * pluginsS;
     179
     180  // set up a CDF midpoint jet definition
     181  #ifdef ENABLE_PLUGIN_CDFCONES
     182  plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
     183  jet_def = fastjet::JetDefinition(plugins);
     184  #else
     185  plugins = NULL;
     186  #endif
     187
     188  // set up a CDF midpoint jet definition
     189  #ifdef ENABLE_PLUGIN_CDFCONES
     190  pluginsS = new fastjet::CDFJetCluPlugin(2,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
     191  jet_defS = fastjet::JetDefinition(pluginsS);
     192  #else
     193  pluginsS = NULL;
     194  #endif
    158195
    159196 
     
    172209      itGen.Reset();
    173210      TrackCentral.clear();
    174       towers.clear();
    175       towersS.clear();
    176211      TSimpleArray<TRootGenParticle> NFCentralQ;
    177 
     212      input_particles.clear();
     213      inclusive_jets.clear();
     214      sorted_jets.clear();
     215      input_particlesS.clear();
     216      inclusive_jetsS.clear();
     217      sorted_jetsS.clear();
    178218
    179219      // Loop over all particles in event
     
    185225          float eta = fabs(particle->Eta);
    186226
    187 if(particle->Status == 1)towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
     227          if(particle->Status == 1)
     228             {
     229               input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     230             }
    188231
    189232          // keeps only final particles, visible by the central detector, including the fiducial volume
     
    224267
    225268                if(pid != pMU && genMomentum.Pt() > DET->PT_TRACKS_MIN)
    226                 towersS.push_back
    227                  (
    228                    PhysicsTower
    229                     (
    230                       LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())
    231                     )
    232                  );
    233                 // all final charged particles
     269                  {
     270                     input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     271                  }
     272               
     273                // all final charged particles
    234274                if(
    235275                     ((rand()%100) < DET->TRACKING_EFF)  &&
     
    251291
    252292      //*****************************
    253       JetCluAlgorithm jetAlgoS(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    254       jetAlgoS.run(towersS, jetsS);
    255 
    256       JetCluAlgorithm jetAlgo(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    257       jetAlgo.run(towers, jets);
    258 
    259 
    260       vector<Cluster>::iterator itJet;
    261       TLorentzVector JETSm(0,0,0,0);
    262       for(itJet = jets.begin(); itJet != jets.end(); ++itJet) {
    263           jetMomentum = itJet->fourVector;
    264           TLorentzVector JET(0,0,0,0);
    265           JET.SetPxPyPzE(jetMomentum.px,jetMomentum.py,jetMomentum.pz,jetMomentum.E);
    266           PairingJet(JETSm,JET,jetsS);
    267           if(JETSm.Pt()>3)
    268             {
    269               elementJet= (RESOLJET*) branchjet->NewEntry();
    270               elementJet->NonSmearePT = JET.Pt();
    271               elementJet->SmearePT = (JETSm.Pt()/JET.Pt());
    272               /*cout<<"valeur obtenue "<<JETSm.Pt()/JET.Pt()<<endl;
    273               cout<<" pt non smeare "<<JET.Pt()<<" phi "<<JET.Phi()<<" eta "<<JET.Eta()<<endl;
    274               cout<<"pt smeare "<<JETSm.Pt()<<" phi "<<JETSm.Phi()<<" eta "<<JETSm.Eta()<<endl;
    275               cout<<"************"<<endl;
     293
     294     double ptmin=1;
     295     if(input_particles.size()!=0)
     296        {
     297          fastjet::ClusterSequence clust_seq(input_particles, jet_def);
     298          inclusive_jets = clust_seq.inclusive_jets(ptmin);
     299          sorted_jets = sorted_by_pt(inclusive_jets);
     300        }
     301
     302     if(input_particlesS.size()!=0)
     303        {
     304          fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
     305          inclusive_jetsS = clust_seqS.inclusive_jets(ptmin);
     306          sorted_jetsS = sorted_by_pt(inclusive_jetsS);
     307        }
     308
     309    TLorentzVector JETSm(0,0,0,0);
     310    for (unsigned int i = 0; i < sorted_jets.size(); i++) {
     311        TLorentzVector JET(0,0,0,0);
     312        JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
     313        PairingJet(JETSm,JET,sorted_jetsS);
     314        if(JETSm.Pt()>1)
     315          {
     316            elementJet= (RESOLJET*) branchjet->NewEntry();
     317            elementJet->NonSmearePT = JET.Pt();
     318            elementJet->SmearePT = (JETSm.Pt()/JET.Pt());
     319            /*cout<<"valeur obtenue "<<JETSm.Pt()/JET.Pt()<<endl;
     320            cout<<" pt non smeare "<<JET.Pt()<<" phi "<<JET.Phi()<<" eta "<<JET.Eta()<<endl;
     321            cout<<"pt smeare "<<JETSm.Pt()<<" phi "<<JETSm.Phi()<<" eta "<<JETSm.Eta()<<endl;
     322            cout<<"************"<<endl;
    276323*/
    277             }
     324          }
    278325
    279326        } // for itJet : loop on all jets
  • trunk/genMakefile.tcl

    r15 r19  
    232232executableDeps {*.cpp}
    233233
    234 dictDeps {DICT} {Utilities/ExRootAnalysis/src/*LinkDef.h}
     234dictDeps {DICT} {Utilities/ExRootAnalysis/src/*LinkDef.h} {src/*LinkDef.h}
    235235
    236236sourceDeps {SOURCE} {src/*.cc} {Utilities/ExRootAnalysis/src/*.cc} {Utilities/Hector/src/*.cc} {Utilities/CDFCones/src/*cc} {Utilities/Fastjet/src/*.cc} {Utilities/Fastjet/plugins/CDFCones/*.cc} {Utilities/Fastjet/plugins/CDFCones/src/*.cc} {Utilities/Fastjet/plugins/SISCone/*.cc} {Utilities/Fastjet/plugins/SISCone/src/*.cc}
Note: See TracChangeset for help on using the changeset viewer.