Changeset 88 in svn for trunk/Resolutions.cpp

Timestamp:

Dec 4, 2008, 3:23:03 PM (16 years ago)

Author:

uid677

Message:

OK for new verson

File:

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

trunk/Resolutions.cpp

@@ -29 +29 @@
 
 #include "interface/SmearUtil.h"
+#include "interface/JetUtils.h"
+#include "interface/BFieldProp.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>
@@ -109 +100 @@
 }
 
-
 int main(int argc, char *argv[])
 {
@@ -133 +123 @@
     return -1;
   }
+
   TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
   outputFile->Close();
@@ -196 +187 @@
   RESOLution *DET = new RESOLution();   
   DET->ReadDataCard(DetDatacard);
-  
-  TLorentzVector genMomentum(0,0,0,0);
+
+  //Jet information
+  JetsUtil *JETRUN = new JetsUtil();
+
+  //Propagation of tracks in the B field
+  TrackPropagation *TRACP = new TrackPropagation();
+   
+  TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
+  TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information
   LorentzVector jetMomentum;
   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;
-
-  vector<fastjet::PseudoJet> input_particlesT;//for FastJet algorithm
-  vector<fastjet::PseudoJet> inclusive_jetsT;
-  vector<fastjet::PseudoJet> sorted_jetsT;
+  vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
+  vector<fastjet::PseudoJet> sorted_jetsGEN;
+  
+  vector<fastjet::PseudoJet> input_particlesReco;//for FastJet algorithm
+  vector<fastjet::PseudoJet> sorted_jetsReco;
 
   vector<PhysicsTower> towers;
-  
-  fastjet::JetDefinition jet_def;
-  fastjet::JetDefinition jet_defS;
-  fastjet::JetDefinition jet_defT;
-  fastjet::JetDefinition::Plugin * plugins;
-  fastjet::JetDefinition::Plugin * pluginsS;
-  fastjet::JetDefinition::Plugin * pluginsT;
-  
-  // set up a CDF midpoint jet definition
-#ifdef ENABLE_PLUGIN_CDFCONES
-  plugins = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->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(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->OVERLAPTHRESHOLD);
-  jet_defS = fastjet::JetDefinition(pluginsS);
-#else
-  pluginsS = NULL;
-#endif
-  
- // set up a CDF midpoint jet definition
- #ifdef ENABLE_PLUGIN_CDFCONES
- pluginsT = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->OVERLAPTHRESHOLD);
- jet_defT = fastjet::JetDefinition(pluginsT);
- #else
- pluginsT = NULL;
- #endif
- 
   
   // Loop over all events
@@ -252 +212 @@
   for(entry = 0; entry < allEntries; ++entry)
     {
-      TLorentzVector PTmisS(0,0,0,0);
-      TLorentzVector PTmis(0,0,0,0);
+      TLorentzVector PTmisReco(0,0,0,0);
+      TLorentzVector PTmisGEN(0,0,0,0);
       treeReader->ReadEntry(entry);
       treeWriter->Clear();
@@ -263 +223 @@
       TSimpleArray<TRootGenParticle> NFCentralQ;
    
-      sorted_jetsS.clear();
-      input_particlesS.clear();
-      inclusive_jetsS.clear();
-
-      sorted_jetsT.clear();
-      input_particlesT.clear();
-      inclusive_jetsT.clear();
-
-      sorted_jets.clear();
-      input_particles.clear();
-      inclusive_jets.clear();
+      input_particlesReco.clear();
+      input_particlesGEN.clear();
       towers.clear();
       
@@ -283 +234 @@
           int pid = abs(particle->PID);
           float eta = fabs(particle->Eta);
-          
+
+          //input generator level particle for jet algorithm      
           if(particle->Status == 1 && eta < DET->MAX_CALO_FWD)
             {
-              input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+              input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
             }
-          
-          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmis = PTmis + genMomentum;
-          // 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) )
-              ) 
-            {
-              //if((pid != pNU1) && (pid != pNU2) && (pid != pNU3))PTmis = PTmis + genMomentum;//ptmis
-
-              Float_t nonS=0;
+
+          //Calculate ETMIS from generated particles
+          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
+
+          if( (particle->Status == 1)    &&
+              ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
+              (fabs(particle->Eta) < DET->MAX_CALO_FWD)
+              )
+            {
+              recoMomentum = genMomentum;          
+              //use of the magnetic field propagation
+              TRACP->Propagation(particle,recoMomentum);
+              float eta=fabs(recoMomentum.Eta());
+
               switch(pid) {
+
               case pE: // all electrons with eta < DET->MAX_CALO_FWD
-                nonS = genMomentum.E();
-                DET->SmearElectron(genMomentum);
+                DET->SmearElectron(recoMomentum);
                 if(eta < DET->MAX_TRACKER){
                 elementElec=(RESOLELEC*) branchelec->NewEntry();
-                elementElec->NonSmeareE = nonS;
-                elementElec->SmeareE = genMomentum.E();}
+                elementElec->E = genMomentum.E();
+                elementElec->SmearedE = recoMomentum.E();}
                 break; // case pE
               case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
-                DET->SmearElectron(genMomentum);
+                DET->SmearElectron(recoMomentum);
                 break; // case pGAMMA
               case pMU: // all muons with eta < DET->MAX_MU
+                DET->SmearMu(recoMomentum);
                 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
-                elementMuon->OneoverNonSmearePT = (1/genMomentum.Pt());
-                DET->SmearMu(genMomentum);
-                elementMuon->OneoverSmearePT = (1/genMomentum.Pt());
+                elementMuon->OverPT = (1/genMomentum.Pt());
+                elementMuon->OverSmearedPT = (1/recoMomentum.Pt());
                 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);
+                DET->SmearHadron(recoMomentum, 0.7);
                 break; // case hadron
               default:   // all other final particles with eta < DET->MAX_CALO_FWD
-                DET->SmearHadron(genMomentum, 1.0);
+                DET->SmearHadron(recoMomentum, 1.0);
                 break;
               } // switch (pid)
-              
-              if(pid != pMU && genMomentum.Et() !=0)
-                {
-                  towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
-                  input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+             
+            //information to reconstruct jets from reconstructed towers
+            int charge=Charge(pid);
+            if(recoMomentum.E() !=0 && pid != pMU) {
+              if(charge == 0 || (charge !=0 && recoMomentum.Pt() >= DET->PT_TRACKS_MIN)){
+                  towers.push_back(PhysicsTower(LorentzVector(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E())));
+                  input_particlesReco.push_back(fastjet::PseudoJet(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E()));
                 }
+            }
               
               // 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
-                 )
-                TrackCentral.push_back(genMomentum);
+                 if(
+                   (genMomentum.E()!=0) &&
+                   (fabs(genMomentum.Eta()) < DET->MAX_TRACKER) &&
+                   (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
+                   ((rand()%100) < DET->TRACKING_EFF)  &&
+                   (charge!=0)
+                   )
+                   {TrackCentral.push_back(genMomentum);}
+
             } // switch
         } // while
-    
+
+     //compute missing transverse energy from calo towers    
      TLorentzVector Att(0.,0.,0.,0.);
      for(unsigned int i=0; i < towers.size(); i++) 
        {
          Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
-         PTmisS = PTmisS + Att;
+         PTmisReco = PTmisReco + Att;
        }
 
       elementEtmis= (ETMIS*) branchetmis->NewEntry();
-      elementEtmis->Et = (PTmis).Pt();
-      elementEtmis->EtSmeare = (PTmisS).Pt()-(PTmis).Pt();
+      elementEtmis->Et = (PTmisGEN).Pt();
+      elementEtmis->SmearedEt = (PTmisReco).Pt()-(PTmisGEN).Pt();
       
       //*****************************
       
-      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);
-        }
-      
+      sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN);
+      sorted_jetsReco=JETRUN->RunJets(input_particlesReco);
+
       TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
 
-      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)
+      TLorentzVector JETreco(0,0,0,0);
+      for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
+        TLorentzVector JETgen(0,0,0,0);
+        JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
+        PairingJet(JETreco,JETgen,sorted_jetsReco);
+        if(JETreco.Pt()>1)
           {
             elementJet= (RESOLJET*) branchjet->NewEntry();
-            elementJet->NonSmearePT = JET.Et();
-            elementJet->SmearePT = JETSm.Et()/JET.Et();
+            elementJet->PT = JETgen.Et();
+            elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
           }
       }
-      for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
+      for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) {
         TLorentzVector JETT(0,0,0,0);
-        JETT.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
+        JETT.SetPxPyPzE(sorted_jetsReco[i].px(),sorted_jetsReco[i].py(),sorted_jetsReco[i].pz(),sorted_jetsReco[i].E());
         if(fabs(JETT.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) 
            {