Fork me on GitHub

Changeset 11 in svn for trunk/Delphes.cpp


Ignore:
Timestamp:
Nov 6, 2008, 3:32:15 PM (16 years ago)
Author:
severine ovyn
Message:

Fastjet added; CDFCones directory has been changed

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Delphes.cpp

    r2 r11  
    1919#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
    2020
    21 #include "Utilities/CDFCones/interface/JetCluAlgorithm.h"
    22 #include "Utilities/CDFCones/interface/MidPointAlgorithm.h"
    23 #include "Utilities/CDFCones/interface/PhysicsTower.h"
    24 #include "Utilities/CDFCones/interface/Cluster.h"
    25 
    2621#include "H_BeamParticle.h"
    2722#include "H_BeamLine.h"
    … …  
    3429
    3530#include "interface/SmearUtil.h"
     31#include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
     32#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
     45
     46#include<vector>
     47#include<iostream>
     48
     49
    3650
    3751using namespace std;
    … …  
    154168  vector<TLorentzVector> TrackCentral; 
    155169  vector<PhysicsTower> towers;
    156   vector<Cluster> jets;
    157 
     170
     171  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
     172  vector<fastjet::PseudoJet> inclusive_jets;
     173  vector<fastjet::PseudoJet> sorted_jets;
     174 
    158175  //Initialisation of Hector
    159176  extern bool relative_energy;
    … …  
    161178  extern int kickers_on;
    162179  kickers_on = 1;         // should always be 1
    163 
    164         // user should provide : (1) optics file for each beamline, and IPname,
    165           // and offset data (s,x) for optical elements
     180 
     181  // user should provide : (1) optics file for each beamline, and IPname,
     182  // and offset data (s,x) for optical elements
    166183  H_BeamLine* beamline1 = new H_BeamLine(1,500.);
    167184  beamline1->fill("data/LHCB1IR5_v6.500.tfs",1,"IP5");
    … …  
    171188  beamline1->add(rp220_1);
    172189  beamline1->add(rp420_1);
    173 
     190 
    174191  H_BeamLine* beamline2 = new H_BeamLine(1,500.);
    175192  beamline2->fill("data/LHCB1IR5_v6.500.tfs",-1,"IP5");
    … …  
    197214      TrackCentral.clear();
    198215      towers.clear();
     216      input_particles.clear();
    199217      TSimpleArray<TRootGenParticle> NFCentralQ;
    200218
    … …  
    266284                      break;
    267285                } // switch (pid)
    268 
     286               
    269287                // all final particles but muons and neutrinos 
    270288                // for calorimetric towers and mission PT
    271289                if(genMomentum.E()!=0) {
    272                         PTmis = PTmis + genMomentum;//ptmis
    273                         if(pid !=pMU) {
    274                               towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
    275                               elementCalo = (TRootCalo*) branchCalo->NewEntry();
    276                               elementCalo->Set(genMomentum);
    277                         }
     290                  PTmis = PTmis + genMomentum;//ptmis
     291                  if(pid !=pMU) {
     292                    towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
     293                    // create a fastjet::PseudoJet with these components and put it onto
     294                    // back of the input_particles vector
     295                    input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     296                    elementCalo = (TRootCalo*) branchCalo->NewEntry();
     297                    elementCalo->Set(genMomentum);
     298                  }
    278299                }
    279 
    280 
     300               
     301               
    281302                // all final charged particles
    282303                if(
    283                      ((rand()%100) < DET->TRACKING_EFF)  &&
    284                      (genMomentum.E()!=0) &&
    285                      (fabs(particle->Eta) < DET->MAX_TRACKER) &&
    286                      (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
    287                      (pid != pGAMMA)    &&
    288                      (pid != pPI0)      &&
    289                      (pid != pK0L)      &&
    290                      (pid != pN)        &&
    291                      (pid != pSIGMA0)   &&
    292                      (pid != pDELTA0)   &&
    293                      (pid != pK0S)   // not charged particles : invisible by tracker
    294                  )
    295                     {
    296                       elementTracks = (TRootTracks*) branchTracks->NewEntry();
    297                       elementTracks->Set(genMomentum);
    298                       TrackCentral.push_back(genMomentum);
    299                     }
    300             } // switch
    301 
     304                   ((rand()%100) < DET->TRACKING_EFF)  &&
     305                   (genMomentum.E()!=0) &&
     306                   (fabs(particle->Eta) < DET->MAX_TRACKER) &&
     307                   (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
     308                   (pid != pGAMMA)    &&
     309                   (pid != pPI0)      &&
     310                   (pid != pK0L)      &&
     311                   (pid != pN)        &&
     312                   (pid != pSIGMA0)   &&
     313                   (pid != pDELTA0)   &&
     314                   (pid != pK0S)   // not charged particles : invisible by tracker
     315                   )
     316                  {
     317                    elementTracks = (TRootTracks*) branchTracks->NewEntry();
     318                    elementTracks->Set(genMomentum);
     319                    TrackCentral.push_back(genMomentum);
     320                  }
     321          } // switch
     322         
    302323          // Forward particles in CASTOR ?
    303 /*        if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
    304                                     && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
    305                
    306 
    307           } // CASTOR
    308 */
     324          /*      if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
     325                  && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
     326                 
     327                 
     328                  } // CASTOR
     329          */
    309330          // Zero degree calorimeter, for forward neutrons and photons
    310331          if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > DET->MIN_ZDC ) {
    311 // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
    312 // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
    313 // !!!!!!!!! pas de smearing ! on garde trop d'info !
    314                 elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
    315                 elementZdc->Set(genMomentum);
    316 
    317                 // time of flight t is t = T + d/[ cos(theta) v ]
    318                 //double tx = acos(particle->Px/particle->Pz);
    319                 //double ty = acos(particle->Py/particle->Pz);
    320                 //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
    321                 //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
    322                 double flight_distance = (DET->ZDC_S - particle->Z*(1E-6));
    323                         // assumes also that the emission angle is so small that 1/(cos theta) = 1
    324                 elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
    325                 elementZdc->side = sign(particle->Eta);
    326 
     332            // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
     333            // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
     334            // !!!!!!!!! pas de smearing ! on garde trop d'info !
     335            elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
     336            elementZdc->Set(genMomentum);
     337           
     338            // time of flight t is t = T + d/[ cos(theta) v ]
     339            //double tx = acos(particle->Px/particle->Pz);
     340            //double ty = acos(particle->Py/particle->Pz);
     341            //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
     342            //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
     343            double flight_distance = (DET->ZDC_S - particle->Z*(1E-6));
     344            // assumes also that the emission angle is so small that 1/(cos theta) = 1
     345            elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
     346            elementZdc->side = sign(particle->Eta);
     347           
    327348          }
    328 
     349         
    329350          // if forward proton
    330351          if( (pid == pP)  && (particle->Status == 1) &&  (fabs(particle->Eta) > DET->MAX_CALO_FWD) )
    331352            {
    332 // !!!!!!!! put here particle->CHARGE and particle->MASS
    333                 H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
    334                 p1.smearAng();
    335                 p1.smearPos();
    336                 p1.setPosition(p1.getX()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,p1.getTY(),0);
    337                 p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
    338 
    339                 H_BeamLine *beamline;
    340                 if(particle->Eta >0) beamline = beamline1;
    341                 else beamline = beamline2;
    342 
    343                 p1.computePath(beamline,1);
    344 
    345                 if(p1.stopped(beamline)) {
    346                     if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
    347                         p1.propagate(DET->RP220_S);
    348                         elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
    349                         elementRP220->X  = (1E-6)*p1.getX();  // [m]
    350                         elementRP220->Y  = (1E-6)*p1.getY();  // [m]
    351                         elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
    352                         elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
    353                         elementRP220->S = p1.getS();          // [m]
    354                         elementRP220->T = -1;           // not yet implemented
    355                         elementRP220->E = p1.getE();    // not yet implemented
    356                         elementRP220->q2 = -1;          // not yet implemented
    357                         elementRP220->side = sign(particle->Eta);
    358 
    359                     } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
    360                         p1.propagate(DET->FP420_S);
    361                         elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
    362                         elementFP420->X  = (1E-6)*p1.getX();  // [m]
    363                         elementFP420->Y  = (1E-6)*p1.getY();  // [m]
    364                         elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
    365                         elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
    366                         elementFP420->S = p1.getS();          // [m]
    367                         elementFP420->T = -1;                 // not yet implemented
    368                         elementFP420->E = p1.getE();          // not yet implemented
    369                         elementFP420->q2 = -1;                // not yet implemented
    370                         elementFP420->side = sign(particle->Eta);
    371                     }
     353              // !!!!!!!! put here particle->CHARGE and particle->MASS
     354              H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
     355              p1.smearAng();
     356              p1.smearPos();
     357              p1.setPosition(p1.getX()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,p1.getTY(),0);
     358              p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
     359             
     360              H_BeamLine *beamline;
     361              if(particle->Eta >0) beamline = beamline1;
     362              else beamline = beamline2;
     363             
     364              p1.computePath(beamline,1);
     365             
     366              if(p1.stopped(beamline)) {
     367                if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
     368                  p1.propagate(DET->RP220_S);
     369                  elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
     370                  elementRP220->X  = (1E-6)*p1.getX();  // [m]
     371                  elementRP220->Y  = (1E-6)*p1.getY();  // [m]
     372                  elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
     373                  elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
     374                  elementRP220->S = p1.getS();          // [m]
     375                  elementRP220->T = -1;                 // not yet implemented
     376                  elementRP220->E = p1.getE();  // not yet implemented
     377                  elementRP220->q2 = -1;                // not yet implemented
     378                  elementRP220->side = sign(particle->Eta);
     379                 
     380                } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
     381                  p1.propagate(DET->FP420_S);
     382                  elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
     383                  elementFP420->X  = (1E-6)*p1.getX();  // [m]
     384                  elementFP420->Y  = (1E-6)*p1.getY();  // [m]
     385                  elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
     386                  elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
     387                  elementFP420->S = p1.getS();        // [m]
     388                  elementFP420->T = -1;                       // not yet implemented
     389                  elementFP420->E = p1.getE();          // not yet implemented
     390                  elementFP420->q2 = -1;                // not yet implemented
     391                  elementFP420->side = sign(particle->Eta);
    372392                }
    373 
    374 //              if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
    375 //                      cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
     393              }
     394             
     395              //                if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
     396              //                        cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
    376397            } // if forward proton
    377 
     398         
    378399        } // while
    379 
     400     
    380401      // computes the Missing Transverse Momentum
    381402      elementEtmis = (TRootETmis*) branchETmis->NewEntry();
    … …  
    384405      elementEtmis->Px = (-PTmis).Px();
    385406      elementEtmis->Py = (-PTmis).Py();
    386 
     407     
    387408      //*****************************
    388       jets.clear();
     409
     410      // we will have four jet definitions, and the first three will be
     411      // plugins
     412      fastjet::JetDefinition jet_def;
     413      fastjet::JetDefinition::Plugin * plugins;
     414     
    389415      switch(DET->JETALGO) {
    390         default:
    391         case 1: {
    392             JetCluAlgorithm jetAlgoC(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    393             jetAlgoC.run(towers, jets);
    394           }
    395           break;
    396 
    397         case 2: {
    398             MidPointAlgorithm jetAlgoM(DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXITERATIONS,DET->M_OVERLAPTHRESHOLD);
    399             jetAlgoM.run(towers, jets);
    400           }
    401           break;
    402 /*
    403         case 3: {
    404                 FastJet
    405           }
    406           break;
    407 */
    408        } // switch
    409 
    410       // Loop on all jets
    411       // Dealing with jets, tau-jets, b-jets
    412       vector<Cluster>::iterator itJet;
    413       for(itJet = jets.begin(); itJet != jets.end(); ++itJet) {
    414           elementJet = (TRootJet*) branchJet->NewEntry();
    415           jetMomentum = itJet->fourVector;
    416           TLorentzVector JET;
    417           JET.SetPxPyPzE(jetMomentum.px,jetMomentum.py,jetMomentum.pz,jetMomentum.E);
    418           elementJet->Set(JET);
    419          
    420           // b-jets
    421           bool btag=false;
    422           if((fabs(JET.Eta()) < DET->MAX_TRACKER && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
    423           elementJet->Btag = btag;
    424 
    425           // Tau jet identification : 1! track and electromagnetic collimation
    426           if(fabs(JET.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) {
    427               double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
    428               if(
    429                   ( Energie_tau_central/JET.E() > DET->TAU_EM_COLLIMATION ) &&
    430                   ( DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 )
    431                 ) {
    432                   elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
    433                   elementTauJet->Set(JET);
    434                } // if tau jet
    435           } // if JET.eta < tracker - tau_cone : Tau jet identification
    436         } // for itJet : loop on all jets
    437 
     416      default:
     417      case 1: {
     418       
     419        // set up a CDF midpoint jet definition
     420        #ifdef ENABLE_PLUGIN_CDFCONES
     421        plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
     422        jet_def = fastjet::JetDefinition(plugins);
     423        #else
     424        plugins = NULL;
     425        #endif
     426       
     427      }
     428        break;
     429       
     430      case 2: {
     431       
     432        // set up a CDF midpoint jet definition
     433        #ifdef ENABLE_PLUGIN_CDFCONES
     434        plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
     435        jet_def = fastjet::JetDefinition(plugins);
     436        #else
     437        plugins = NULL;
     438        #endif
     439      }
     440        break;
     441       
     442      case 3: {
     443        // set up a siscone jet definition
     444        #ifdef ENABLE_PLUGIN_SISCONE
     445        int npass = 0;               // do infinite number of passes
     446        double protojet_ptmin = 0.0; // use all protojets
     447        plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
     448        jet_def = fastjet::JetDefinition(plugins);
     449        #else
     450        plugins = NULL;
     451        #endif
     452      }
     453        break;
     454       
     455      case 4: {
     456        jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
     457        //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
     458        //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
     459      }
     460        break;
     461      }
     462      // run the jet clustering with the above jet definition
     463      if(input_particles.size()!=0)
     464        {
     465          fastjet::ClusterSequence clust_seq(input_particles, jet_def);
     466         
     467         
     468          // extract the inclusive jets with pt > 5 GeV
     469          double ptmin = 5.0;
     470          inclusive_jets = clust_seq.inclusive_jets(ptmin);
     471         
     472          // sort jets into increasing pt
     473          sorted_by_pt(inclusive_jets);
     474        }
     475     
     476      for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
     477        elementJet = (TRootJet*) branchJet->NewEntry();
     478        TLorentzVector JET;
     479        JET.SetPxPyPzE(inclusive_jets[i].px(),inclusive_jets[i].py(),inclusive_jets[i].pz(),inclusive_jets[i].E());
     480        elementJet->Set(JET);
     481        // b-jets
     482        bool btag=false;
     483        if((fabs(JET.Eta()) < DET->MAX_TRACKER && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
     484        elementJet->Btag = btag;
     485       
     486        // Tau jet identification : 1! track and electromagnetic collimation
     487        if(fabs(JET.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) {
     488          double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
     489          if(
     490             ( Energie_tau_central/JET.E() > DET->TAU_EM_COLLIMATION ) &&
     491             ( DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 )
     492             ) {
     493            elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
     494            elementTauJet->Set(JET);
     495          } // if tau jet
     496        } // if JET.eta < tracker - tau_cone : Tau jet identification
     497      } // for itJet : loop on all jets
     498     
    438499      treeWriter->Fill();
    439500      // Add here the trigger
Note: See TracChangeset for help on using the changeset viewer.