Changeset 19 in svn for trunk/Resolutions.cpp

Timestamp:

Nov 6, 2008, 5:48:51 PM (16 years ago)

Author:

severine ovyn

Message:

plugins at the begining

File:

• trunk/Resolutions.cpp (modified) (8 diffs)

trunk/Resolutions.cpp

@@ -19 +19 @@
 #include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
 
-#include "Utilities/CDFCones/interface/JetCluAlgorithm.h"
-#include "Utilities/CDFCones/interface/MidPointAlgorithm.h"
-#include "Utilities/CDFCones/interface/PhysicsTower.h"
-#include "Utilities/CDFCones/interface/Cluster.h"
+#include "H_BeamParticle.h"
+#include "H_BeamLine.h"
+#include "H_RomanPot.h"
 
 #include "interface/DataConverter.h"
@@ -28 +27 @@
 #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>
+
 #include "interface/TreeClasses.h"
-
 using namespace std;
 
@@ -37 +54 @@
 
 
-void PairingJet(TLorentzVector &JETSm, TLorentzVector JET, vector<Cluster> jetsS)
+void PairingJet(TLorentzVector &JETSm, TLorentzVector JET,  vector<fastjet::PseudoJet> sorted_jetsS)
 {
   JETSm.SetPxPyPzE(0,0,0,0);
-  vector<Cluster>::iterator itJetS;
   float deltaRtest=5000;
-  LorentzVector jetMomentumS;
-  for(itJetS = jetsS.begin(); itJetS != jetsS.end(); ++itJetS) 
-     {
-       jetMomentumS = itJetS->fourVector;
-       TLorentzVector Att;
-       Att.SetPxPyPzE(jetMomentumS.px,jetMomentumS.py,jetMomentumS.pz,jetMomentumS.E);
-       if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
-            {
-              deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
-              if(deltaRtest < 0.25)
-                {
-                  JETSm = Att;
-                }
-            }
-     }
+  for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
+      TLorentzVector Att;
+      Att.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
+      if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
+        {
+          deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
+          if(deltaRtest < 0.25)
+             {
+               JETSm = Att;
+             }
+        }
+    }
 }
 
@@ -148 +161 @@
   DET->ReadDataCard(DetDatacard);
   
-  
   TLorentzVector genMomentum(0,0,0,0);
   LorentzVector jetMomentum;
-  vector<TLorentzVector> TrackCentral;  
-  vector<PhysicsTower> towersS;
-  vector<Cluster> jetsS; 
-  vector<PhysicsTower> towers;
-  vector<Cluster> jets;
+  vector<TLorentzVector> TrackCentral;
+
+  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
+  vector<fastjet::PseudoJet> inclusive_jets;
+  vector<fastjet::PseudoJet> sorted_jets;
+
+  vector<fastjet::PseudoJet> input_particlesS;//for FastJet algorithm
+  vector<fastjet::PseudoJet> inclusive_jetsS;
+  vector<fastjet::PseudoJet> sorted_jetsS;
+
+  fastjet::JetDefinition jet_def;
+  fastjet::JetDefinition jet_defS;
+  fastjet::JetDefinition::Plugin * plugins;
+  fastjet::JetDefinition::Plugin * pluginsS;
+
+  // 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
+
+  // set up a CDF midpoint jet definition
+  #ifdef ENABLE_PLUGIN_CDFCONES
+  pluginsS = new fastjet::CDFJetCluPlugin(2,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
+  jet_defS = fastjet::JetDefinition(pluginsS);
+  #else
+  pluginsS = NULL;
+  #endif
 
  
@@ -172 +209 @@
       itGen.Reset();
       TrackCentral.clear();
-      towers.clear();
-      towersS.clear();
       TSimpleArray<TRootGenParticle> NFCentralQ;
-
+      input_particles.clear();
+      inclusive_jets.clear();
+      sorted_jets.clear();
+      input_particlesS.clear();
+      inclusive_jetsS.clear();
+      sorted_jetsS.clear();
 
       // Loop over all particles in event
@@ -185 +225 @@
           float eta = fabs(particle->Eta);
 
-if(particle->Status == 1)towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
+          if(particle->Status == 1)
+             {
+               input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+             }
 
           // keeps only final particles, visible by the central detector, including the fiducial volume
@@ -224 +267 @@
 
                 if(pid != pMU && genMomentum.Pt() > DET->PT_TRACKS_MIN)
-                towersS.push_back
-                 (
-                   PhysicsTower
-                    (
-                      LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())
-                    )
-                 );
-                // all final charged particles
+                  {
+                     input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+                  }
+                
+                // all final charged particles
                 if(
                      ((rand()%100) < DET->TRACKING_EFF)  &&
@@ -251 +291 @@
 
       //*****************************
-      JetCluAlgorithm jetAlgoS(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
-      jetAlgoS.run(towersS, jetsS);
-
-      JetCluAlgorithm jetAlgo(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
-      jetAlgo.run(towers, jets);
-
-
-      vector<Cluster>::iterator itJet;
-      TLorentzVector JETSm(0,0,0,0);
-      for(itJet = jets.begin(); itJet != jets.end(); ++itJet) {
-          jetMomentum = itJet->fourVector;
-          TLorentzVector JET(0,0,0,0);
-          JET.SetPxPyPzE(jetMomentum.px,jetMomentum.py,jetMomentum.pz,jetMomentum.E);
-          PairingJet(JETSm,JET,jetsS);
-          if(JETSm.Pt()>3)
-            {
-              elementJet= (RESOLJET*) branchjet->NewEntry();
-              elementJet->NonSmearePT = JET.Pt();
-              elementJet->SmearePT = (JETSm.Pt()/JET.Pt());
-              /*cout<<"valeur obtenue "<<JETSm.Pt()/JET.Pt()<<endl;
-              cout<<" pt non smeare "<<JET.Pt()<<" phi "<<JET.Phi()<<" eta "<<JET.Eta()<<endl;
-              cout<<"pt smeare "<<JETSm.Pt()<<" phi "<<JETSm.Phi()<<" eta "<<JETSm.Eta()<<endl;
-              cout<<"************"<<endl;
+
+     double ptmin=1;
+     if(input_particles.size()!=0)
+        {
+          fastjet::ClusterSequence clust_seq(input_particles, jet_def);
+          inclusive_jets = clust_seq.inclusive_jets(ptmin);
+          sorted_jets = sorted_by_pt(inclusive_jets);
+        }
+
+     if(input_particlesS.size()!=0)
+        {
+          fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
+          inclusive_jetsS = clust_seqS.inclusive_jets(ptmin);
+          sorted_jetsS = sorted_by_pt(inclusive_jetsS);
+        }
+
+    TLorentzVector JETSm(0,0,0,0);
+    for (unsigned int i = 0; i < sorted_jets.size(); i++) {
+        TLorentzVector JET(0,0,0,0);
+        JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
+        PairingJet(JETSm,JET,sorted_jetsS);
+        if(JETSm.Pt()>1)
+          {
+            elementJet= (RESOLJET*) branchjet->NewEntry();
+            elementJet->NonSmearePT = JET.Pt();
+            elementJet->SmearePT = (JETSm.Pt()/JET.Pt());
+            /*cout<<"valeur obtenue "<<JETSm.Pt()/JET.Pt()<<endl;
+            cout<<" pt non smeare "<<JET.Pt()<<" phi "<<JET.Phi()<<" eta "<<JET.Eta()<<endl;
+            cout<<"pt smeare "<<JETSm.Pt()<<" phi "<<JETSm.Phi()<<" eta "<<JETSm.Eta()<<endl;
+            cout<<"************"<<endl;
 */
-            }
+          }
 
         } // for itJet : loop on all jets