Fork me on GitHub

Changeset 88 in svn


Ignore:
Timestamp:
Dec 4, 2008, 3:23:03 PM (16 years ago)
Author:
uid677
Message:

OK for new verson

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Resolutions.cpp

    r43 r88  
    2929
    3030#include "interface/SmearUtil.h"
     31#include "interface/JetUtils.h"
     32#include "interface/BFieldProp.h"
     33
    3134#include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
    3235#include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh"
    33 
    34 // get info on how fastjet was configured
    35 #include "Utilities/Fastjet/include/fastjet/config.h"
    36 
    37 // make sure we have what is needed
    38 #ifdef ENABLE_PLUGIN_SISCONE
    39 #  include "Utilities/Fastjet/plugins/SISCone/SISConePlugin.hh"
    40 #endif
    41 #ifdef ENABLE_PLUGIN_CDFCONES
    42 #  include "Utilities/Fastjet/plugins/CDFCones/CDFMidPointPlugin.hh"
    43 #  include "Utilities/Fastjet/plugins/CDFCones/CDFJetCluPlugin.hh"
    44 #endif
    4536
    4637#include<vector>
     
    109100}
    110101
    111 
    112102int main(int argc, char *argv[])
    113103{
     
    133123    return -1;
    134124  }
     125
    135126  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
    136127  outputFile->Close();
     
    196187  RESOLution *DET = new RESOLution();   
    197188  DET->ReadDataCard(DetDatacard);
    198  
    199   TLorentzVector genMomentum(0,0,0,0);
     189
     190  //Jet information
     191  JetsUtil *JETRUN = new JetsUtil();
     192
     193  //Propagation of tracks in the B field
     194  TrackPropagation *TRACP = new TrackPropagation();
     195   
     196  TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
     197  TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information
    200198  LorentzVector jetMomentum;
    201199  vector<TLorentzVector> TrackCentral;
    202200 
    203   vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
    204   vector<fastjet::PseudoJet> inclusive_jets;
    205   vector<fastjet::PseudoJet> sorted_jets;
    206  
    207   vector<fastjet::PseudoJet> input_particlesS;//for FastJet algorithm
    208   vector<fastjet::PseudoJet> inclusive_jetsS;
    209   vector<fastjet::PseudoJet> sorted_jetsS;
    210 
    211   vector<fastjet::PseudoJet> input_particlesT;//for FastJet algorithm
    212   vector<fastjet::PseudoJet> inclusive_jetsT;
    213   vector<fastjet::PseudoJet> sorted_jetsT;
     201  vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
     202  vector<fastjet::PseudoJet> sorted_jetsGEN;
     203 
     204  vector<fastjet::PseudoJet> input_particlesReco;//for FastJet algorithm
     205  vector<fastjet::PseudoJet> sorted_jetsReco;
    214206
    215207  vector<PhysicsTower> towers;
    216  
    217   fastjet::JetDefinition jet_def;
    218   fastjet::JetDefinition jet_defS;
    219   fastjet::JetDefinition jet_defT;
    220   fastjet::JetDefinition::Plugin * plugins;
    221   fastjet::JetDefinition::Plugin * pluginsS;
    222   fastjet::JetDefinition::Plugin * pluginsT;
    223  
    224   // set up a CDF midpoint jet definition
    225 #ifdef ENABLE_PLUGIN_CDFCONES
    226   plugins = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->OVERLAPTHRESHOLD);
    227   jet_def = fastjet::JetDefinition(plugins);
    228 #else
    229   plugins = NULL;
    230 #endif
    231  
    232   // set up a CDF midpoint jet definition
    233 #ifdef ENABLE_PLUGIN_CDFCONES
    234   pluginsS = new fastjet::CDFJetCluPlugin(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->OVERLAPTHRESHOLD);
    235   jet_defS = fastjet::JetDefinition(pluginsS);
    236 #else
    237   pluginsS = NULL;
    238 #endif
    239  
    240  // set up a CDF midpoint jet definition
    241  #ifdef ENABLE_PLUGIN_CDFCONES
    242  pluginsT = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->OVERLAPTHRESHOLD);
    243  jet_defT = fastjet::JetDefinition(pluginsT);
    244  #else
    245  pluginsT = NULL;
    246  #endif
    247  
    248208 
    249209  // Loop over all events
     
    252212  for(entry = 0; entry < allEntries; ++entry)
    253213    {
    254       TLorentzVector PTmisS(0,0,0,0);
    255       TLorentzVector PTmis(0,0,0,0);
     214      TLorentzVector PTmisReco(0,0,0,0);
     215      TLorentzVector PTmisGEN(0,0,0,0);
    256216      treeReader->ReadEntry(entry);
    257217      treeWriter->Clear();
     
    263223      TSimpleArray<TRootGenParticle> NFCentralQ;
    264224   
    265       sorted_jetsS.clear();
    266       input_particlesS.clear();
    267       inclusive_jetsS.clear();
    268 
    269       sorted_jetsT.clear();
    270       input_particlesT.clear();
    271       inclusive_jetsT.clear();
    272 
    273       sorted_jets.clear();
    274       input_particles.clear();
    275       inclusive_jets.clear();
     225      input_particlesReco.clear();
     226      input_particlesGEN.clear();
    276227      towers.clear();
    277228     
     
    283234          int pid = abs(particle->PID);
    284235          float eta = fabs(particle->Eta);
    285          
     236
     237          //input generator level particle for jet algorithm     
    286238          if(particle->Status == 1 && eta < DET->MAX_CALO_FWD)
    287239            {
    288               input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     240              input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
    289241            }
    290          
    291           if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmis = PTmis + genMomentum;
    292           // keeps only final particles, visible by the central detector, including the fiducial volume
    293           // the ordering of conditions have been optimised for speed : put first the STATUS condition
    294           if( (particle->Status == 1)    &&
    295               (
    296                (pid == pMU  && eta < DET->MAX_MU) ||
    297                (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) )
    298               )
    299             {
    300               //if((pid != pNU1) && (pid != pNU2) && (pid != pNU3))PTmis = PTmis + genMomentum;//ptmis
    301 
    302               Float_t nonS=0;
     242
     243          //Calculate ETMIS from generated particles
     244          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
     245
     246          if( (particle->Status == 1)    &&
     247              ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
     248              (fabs(particle->Eta) < DET->MAX_CALO_FWD)
     249              )
     250            {
     251              recoMomentum = genMomentum;         
     252              //use of the magnetic field propagation
     253              TRACP->Propagation(particle,recoMomentum);
     254              float eta=fabs(recoMomentum.Eta());
     255
    303256              switch(pid) {
     257
    304258              case pE: // all electrons with eta < DET->MAX_CALO_FWD
    305                 nonS = genMomentum.E();
    306                 DET->SmearElectron(genMomentum);
     259                DET->SmearElectron(recoMomentum);
    307260                if(eta < DET->MAX_TRACKER){
    308261                elementElec=(RESOLELEC*) branchelec->NewEntry();
    309                 elementElec->NonSmeareE = nonS;
    310                 elementElec->SmeareE = genMomentum.E();}
     262                elementElec->E = genMomentum.E();
     263                elementElec->SmearedE = recoMomentum.E();}
    311264                break; // case pE
    312265              case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
    313                 DET->SmearElectron(genMomentum);
     266                DET->SmearElectron(recoMomentum);
    314267                break; // case pGAMMA
    315268              case pMU: // all muons with eta < DET->MAX_MU
     269                DET->SmearMu(recoMomentum);
    316270                elementMuon = (RESOLMUON*) branchmuon->NewEntry();
    317                 elementMuon->OneoverNonSmearePT = (1/genMomentum.Pt());
    318                 DET->SmearMu(genMomentum);
    319                 elementMuon->OneoverSmearePT = (1/genMomentum.Pt());
     271                elementMuon->OverPT = (1/genMomentum.Pt());
     272                elementMuon->OverSmearedPT = (1/recoMomentum.Pt());
    320273                break; // case pMU
    321274              case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
    322275              case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
    323                 DET->SmearHadron(genMomentum, 0.7);
     276                DET->SmearHadron(recoMomentum, 0.7);
    324277                break; // case hadron
    325278              default:   // all other final particles with eta < DET->MAX_CALO_FWD
    326                 DET->SmearHadron(genMomentum, 1.0);
     279                DET->SmearHadron(recoMomentum, 1.0);
    327280                break;
    328281              } // switch (pid)
    329              
    330               if(pid != pMU && genMomentum.Et() !=0)
    331                 {
    332                   towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
    333                   input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     282             
     283            //information to reconstruct jets from reconstructed towers
     284            int charge=Charge(pid);
     285            if(recoMomentum.E() !=0 && pid != pMU) {
     286              if(charge == 0 || (charge !=0 && recoMomentum.Pt() >= DET->PT_TRACKS_MIN)){
     287                  towers.push_back(PhysicsTower(LorentzVector(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E())));
     288                  input_particlesReco.push_back(fastjet::PseudoJet(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E()));
    334289                }
     290            }
    335291             
    336292              // all final charged particles
    337               if(
    338                  ((rand()%100) < DET->TRACKING_EFF)  &&
    339                  (genMomentum.E()!=0) &&
    340                  (fabs(particle->Eta) < DET->MAX_TRACKER) &&
    341                  (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
    342                  (pid != pGAMMA)    &&
    343                  (pid != pPI0)      &&
    344                  (pid != pK0L)      &&
    345                  (pid != pN)        &&
    346                  (pid != pSIGMA0)   &&
    347                  (pid != pDELTA0)   &&
    348                  (pid != pK0S)   // not charged particles : invisible by tracker
    349                  )
    350                 TrackCentral.push_back(genMomentum);
     293                 if(
     294                   (genMomentum.E()!=0) &&
     295                   (fabs(genMomentum.Eta()) < DET->MAX_TRACKER) &&
     296                   (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
     297                   ((rand()%100) < DET->TRACKING_EFF)  &&
     298                   (charge!=0)
     299                   )
     300                   {TrackCentral.push_back(genMomentum);}
     301
    351302            } // switch
    352303        } // while
    353    
     304
     305     //compute missing transverse energy from calo towers   
    354306     TLorentzVector Att(0.,0.,0.,0.);
    355307     for(unsigned int i=0; i < towers.size(); i++)
    356308       {
    357309         Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
    358          PTmisS = PTmisS + Att;
     310         PTmisReco = PTmisReco + Att;
    359311       }
    360312
    361313      elementEtmis= (ETMIS*) branchetmis->NewEntry();
    362       elementEtmis->Et = (PTmis).Pt();
    363       elementEtmis->EtSmeare = (PTmisS).Pt()-(PTmis).Pt();
     314      elementEtmis->Et = (PTmisGEN).Pt();
     315      elementEtmis->SmearedEt = (PTmisReco).Pt()-(PTmisGEN).Pt();
    364316     
    365317      //*****************************
    366318     
    367       double ptmin=1;
    368       if(input_particles.size()!=0)
    369         {
    370           fastjet::ClusterSequence clust_seq(input_particles, jet_def);
    371           inclusive_jets = clust_seq.inclusive_jets(ptmin);
    372           sorted_jets = sorted_by_pt(inclusive_jets);
    373         }
    374      
    375       if(input_particlesS.size()!=0)
    376         {
    377           fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
    378           inclusive_jetsS = clust_seqS.inclusive_jets(ptmin);
    379           sorted_jetsS = sorted_by_pt(inclusive_jetsS);
    380         }
    381      
     319      sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN);
     320      sorted_jetsReco=JETRUN->RunJets(input_particlesReco);
     321
    382322      TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
    383323
    384       TLorentzVector JETSm(0,0,0,0);
    385       for (unsigned int i = 0; i < sorted_jets.size(); i++) {
    386         TLorentzVector JET(0,0,0,0);
    387         JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
    388         PairingJet(JETSm,JET,sorted_jetsS);
    389         if(JETSm.Pt()>1)
     324      TLorentzVector JETreco(0,0,0,0);
     325      for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
     326        TLorentzVector JETgen(0,0,0,0);
     327        JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
     328        PairingJet(JETreco,JETgen,sorted_jetsReco);
     329        if(JETreco.Pt()>1)
    390330          {
    391331            elementJet= (RESOLJET*) branchjet->NewEntry();
    392             elementJet->NonSmearePT = JET.Et();
    393             elementJet->SmearePT = JETSm.Et()/JET.Et();
     332            elementJet->PT = JETgen.Et();
     333            elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
    394334          }
    395335      }
    396       for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
     336      for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) {
    397337        TLorentzVector JETT(0,0,0,0);
    398         JETT.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
     338        JETT.SetPxPyPzE(sorted_jetsReco[i].px(),sorted_jetsReco[i].py(),sorted_jetsReco[i].pz(),sorted_jetsReco[i].E());
    399339        if(fabs(JETT.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS))
    400340           {
  • trunk/interface/TreeClasses.h

    r17 r88  
    1313{
    1414public:
    15 Float_t SmearePT;
    16 Float_t NonSmearePT;
     15Float_t SmearedPT;
     16Float_t PT;
    1717
    1818  static TCompare *fgCompare; //!
     
    3030{
    3131public:
    32 Float_t SmeareE;
    33 Float_t NonSmeareE;
     32Float_t E;
     33Float_t SmearedE;
    3434
    3535  static TCompare *fgCompare; //!
     
    4646{
    4747public:
    48 Float_t OneoverSmearePT;
    49 Float_t OneoverNonSmearePT;
     48Float_t OverSmearedPT;
     49Float_t OverPT;
    5050
    5151  static TCompare *fgCompare; //!
     
    7777public:
    7878Float_t Et;
    79 Float_t EtSmeare;
     79Float_t SmearedEt;
    8080
    8181 static TCompare *fgCompare; //!
Note: See TracChangeset for help on using the changeset viewer.