Changeset 447 in svn

Timestamp:

Jun 19, 2009, 7:07:04 PM (15 years ago)

Author:

Xavier Rouby

Message:

update resolutions CMS/ATLAS

Location:

trunk

Files:

• Resolutions.cpp (modified) (14 diffs)
• Resolutions_ATLAS.cpp (added)

trunk/Resolutions.cpp

@@ -33 +33 @@
 
 /// \file Resolution.cpp
-/// \brief Prepares the resolution calculation
+/// \brief Resolution for CMS
 
 #include "TChain.h"
@@ -111 +111 @@
 
 
-void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet, const float dR=0.25)
+void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet)
 {
   JETSm.SetPxPyPzE(0,0,0,0);
@@ -125 +125 @@
         {
           deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
-          if(deltaRtest < dR)
+          if(deltaRtest < 0.25)
              {
                JETSm = Att;
@@ -203 +203 @@
   if(argc != 3) {
       cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
-      cout << " input_file - input file in Delphes-root format," << endl;
+      cout << " input_list - list of files in root format," << endl;
       cout << " output_file - output file." << endl;
       exit(1);
@@ -220 +220 @@
 
 
-  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE");// Creates the file, but should be closed just after
+  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
   outputFile->Close();
   
@@ -257 +257 @@
   
   RESOLution *DET = new RESOLution();
-  /*
-  string detectorcard = "data/DetectorCard_CMS.dat";
-  const float dR_jetpairing = 0.25;
-  const float jet_pt_cut = 1;
-  */
-  string detectorcard = "data/DetectorCard_ATLAS.dat";
-  const float dR_jetpairing = 0.2;
-  const float jet_pt_cut = 7;
-  DET->ReadDataCard(detectorcard);
+  DET->ReadDataCard("data/DetectorCard_CMS.dat");
 
 
   //Jet information
-  JetsUtil *JETRUN = new JetsUtil(detectorcard);
+  JetsUtil *JETRUN = new JetsUtil("data/DetectorCard_CMS.dat");
 
   TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
@@ -277 +269 @@
   vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
   vector<fastjet::PseudoJet> sorted_jetsGEN;
+ 
   vector<int> NTrackJet;
+ 
   vector<TLorentzVector> towers;
  
@@ -315 +309 @@
           //Calculate ETMIS from generated particles
           if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
-
-          //Electrons and muons  
+         
           if( (particle->Status == 1)    &&
             ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
@@ -335 +328 @@
                  case pMU: // all muons with eta < DET->MAX_MU
                     PairingMuon(recoMomentum,genMomentum,branchMuon);
-                    if(recoMomentum.E()!=0){
+                    if(recoMomentum.E() !=0){
                      elementMuon = (RESOLMUON*) branchmuon->NewEntry();
-                     elementMuon->OverPT = 1./genMomentum.Pt();
-                     elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
+                     elementMuon->OverPT = (1/genMomentum.Pt());
+                     elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
                    break; // case pMU
                  default:
@@ -369 +362 @@
       elementEtmis->Ex = (-PTmisGEN).Px();
       elementEtmis->SEt = ScalarEt;
+
       elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
       elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
@@ -376 +370 @@
 
       TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
+
       TLorentzVector JETreco(0,0,0,0);
       for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
@@ -381 +376 @@
         JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
         PairingJet(JETreco,JETgen,branchJet);
-        if(JETreco.Pt()>jet_pt_cut)
+        if(JETreco.Pt()>1)
           {
             elementJet= (RESOLJET*) branchjet->NewEntry();
             elementJet->PT = JETgen.Et();
             elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
-            elementJet->E = JETgen.E();
-            elementJet->dE = (JETreco.E()-JETgen.E())/JETgen.E() ;
-            elementJet->dE2 = pow( (JETreco.E()-JETgen.E())/JETgen.E() , 2.) ;
           }
       }
@@ -407 +399 @@
                   {
                     elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
-                    elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
+                    elementTaujet->EnergieCen = (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E());
                     elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
                     if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95 
@@ -417 +409 @@
         
       } // for itJet : loop on all jets
+//cout<<"i"<<endl;
       
       treeWriter->Fill();
     } // Loop over all events
   treeWriter->Write();
-
-  cout << detectorcard << " has been used.\n";
+float frac = numTauRec/numTau;  
+cout<<numTauRec<<endl;
+cout<<numTau<<endl;
+
   cout << "** Exiting..." << endl;
+  cout<<frac<<endl;
+
   
   delete treeWriter;