Changeset 55 in svn for trunk/Delphes.cpp

Timestamp:

Nov 27, 2008, 3:05:11 PM (16 years ago)

Author:

severine ovyn

Message:

1file ->several

File:

• trunk/Delphes.cpp (modified) (13 diffs)

trunk/Delphes.cpp

@@ -19 +19 @@
 #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"
@@ -29 +25 @@
 
 #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/BFieldProp.h"
+#include "interface/TriggerUtil.h"
+#include "interface/VeryForward.h"
+#include "interface/JetUtils.h"
+
+#include <vector>
+#include <iostream>
 
 using namespace std;
@@ -68 +52 @@
 {
   int appargc = 2;
-  char *appName = "JetsSim";
+  char *appName = "Delphes";
   char *appargv[] = {appName, "-b"};
   TApplication app(appName, &appargc, appargv);
@@ -103 +87 @@
   string DetDatacard("");
   if(argc==4)  DetDatacard =argv[3];
+
+  //Smearing information
   RESOLution *DET = new RESOLution();
   DET->ReadDataCard(DetDatacard);
   DET->Logfile(LogName);
 
-  todo(LogName.c_str());
+  //Trigger information
+  Trigger *TRIG = new Trigger();
+  TRIG->TriggerReader("data/trigger.dat");
+
+  //Propagation of tracks in the B field
+  TrackPropagation *TRACP = new TrackPropagation(); 
+
+  //Jet information
+  JetsUtil *JETRUN = new JetsUtil();
+
+  //VFD information
+  VeryForward * VFD = new VeryForward();
+
+  //todo(LogName.c_str());
  
   DataConverter *converter=0;
@@ -163 +162 @@
   TRootMuon *elementMu;
   TRootPhoton *elementPhoton;
-  TRootJet *elementJet;
-  TRootTauJet *elementTauJet;
   TRootTracks *elementTracks;
   TRootCalo *elementCalo;
-  TRootZdcHits *elementZdc;
-  TRootRomanPotHits *elementRP220, *elementFP420;
   
   TLorentzVector genMomentum(0,0,0,0);
   TLorentzVector genMomentumCalo(0,0,0,0);
   LorentzVector jetMomentum;
+ 
+  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
+  vector<fastjet::PseudoJet> sorted_jets;
+
   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;
@@ -185 +181 @@
   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->SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->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->SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXITERATIONS,DET->OVERLAPTHRESHOLD);
-    jet_def = fastjet::JetDefinition(plugins);
-#else
-    plugins = NULL;
-#endif
-  }
-    break;
-  case 3: {
-    // set up a siscone jet definition
-#ifdef ENABLE_PLUGIN_SISCONE
-    plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->OVERLAPTHRESHOLD,DET->NPASS, DET->PROTOJET_PTMIN);
-    jet_def = fastjet::JetDefinition(plugins);
-#else
-    plugins = NULL;
-#endif
-  }
-    break;
-    
-  case 4: {
-    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
-  }
-    break;
-    
-  case 5: {
-    jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,DET->CONERADIUS);
-  }
-    break;
-    
-  case 6: {
-    jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,DET->CONERADIUS);
-  }
-    break;
-  }
-  
+
+
   // Loop over all events
   Long64_t entry, allEntries = treeReader->GetEntries();
@@ -286 +203 @@
       towers.clear();
       input_particles.clear();
-      inclusive_jets.clear();
-      sorted_jets.clear();
       
       // Loop over all particles in event
@@ -294 +209 @@
           genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
           
-          int pid = abs(particle->PID);
-          float eta = fabs(particle->Eta);
-          
+          int pid = abs(particle->PID);
+          float eta = fabs(genMomentum.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)
@@ -312 +227 @@
           // 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)    &&
               ( 
@@ -318 +235 @@
                 )
               ) { 
-            switch(pid) {
+          //  TRACP->Propagation(particle,genMomentum);
+            eta = fabs(genMomentum.Eta());
+            switch(pid) {
               
             case pE: // all electrons with eta < DET->MAX_CALO_FWD
@@ -363 +282 @@
             }
             
-            
             // all final charged particles
             if(
                ((rand()%100) < DET->TRACKING_EFF)  &&
                (genMomentum.E()!=0) &&
-               (fabs(particle->Eta) < DET->MAX_TRACKER) && 
+               (fabs(genMomentum.Eta()) < DET->MAX_TRACKER) && 
                (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
                (pid != pGAMMA)    && 
@@ -386 +304 @@
           } // 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
+          VFD->ZDC(treeWriter,branchZDC,particle);
+          VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
           
         } // while
@@ -497 +342 @@
       //*****************************
       
-      // 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++) {
-        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 = (TRootJet*) branchJet->NewEntry();
-            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
-      
+      sorted_jets=JETRUN->RunJets(input_particles);
+      JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ);
+      JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral);
+
+
       treeWriter->Fill();
+
       // Add here the trigger
       // Should test all the trigger table on the event, based on reconstructed objects
+      
+ //     Trigger *TRIG = new Trigger();
+ //     TRIG->TriggerReader("data/trigger.dat");
+     
     } // Loop over all events
+
   treeWriter->Write();