Changeset 445 in svn

Timestamp:

Jun 19, 2009, 6:59:40 PM (15 years ago)

Author:

Xavier Rouby

Message:

modifs pour ATLAS resolution

File:

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

trunk/Resolutions.cpp

@@ -32 +32 @@
 ***********************************************************************/
 
-/// \file Smearing.cpp
-/// \brief executable for the FastSim
+/// \file Resolution.cpp
+/// \brief Prepares the resolution calculation
 
 #include "TChain.h"
@@ -111 +111 @@
 
 
-void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet)
+void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet, const float dR=0.25)
 {
   JETSm.SetPxPyPzE(0,0,0,0);
@@ -125 +125 @@
         {
           deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
-          if(deltaRtest < 0.25)
+          if(deltaRtest < dR)
              {
                JETSm = Att;
@@ -203 +203 @@
   if(argc != 3) {
       cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
-      cout << " input_list - list of files in root format," << endl;
+      cout << " input_file - input file in Delphes-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();
-  DET->ReadDataCard("data/DetectorCard_CMS.dat");
+  /*
+  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);
 
 
   //Jet information
-  JetsUtil *JETRUN = new JetsUtil("data/DetectorCard_CMS.dat");
+  JetsUtil *JETRUN = new JetsUtil(detectorcard);
 
   TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
@@ -269 +277 @@
   vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
   vector<fastjet::PseudoJet> sorted_jetsGEN;
- 
   vector<int> NTrackJet;
- 
   vector<TLorentzVector> towers;
  
@@ -309 +315 @@
           //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)) &&
@@ -328 +335 @@
                  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:
@@ -362 +369 @@
       elementEtmis->Ex = (-PTmisGEN).Px();
       elementEtmis->SEt = ScalarEt;
-
       elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
       elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
@@ -370 +376 @@
 
       TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
-
       TLorentzVector JETreco(0,0,0,0);
       for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
@@ -376 +381 @@
         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()>1)
+        if(JETreco.Pt()>jet_pt_cut)
           {
             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.) ;
           }
       }
@@ -399 +407 @@
                   {
                     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 
@@ -409 +417 @@
         
       } // for itJet : loop on all jets
-//cout<<"i"<<endl;
       
       treeWriter->Fill();
     } // Loop over all events
   treeWriter->Write();
-float frac = numTauRec/numTau;  
-cout<<numTauRec<<endl;
-cout<<numTau<<endl;
-
+
+  cout << detectorcard << " has been used.\n";
   cout << "** Exiting..." << endl;
-  cout<<frac<<endl;
-
   
   delete treeWriter;