Utils/src/PhotonTools.cc

// $Id: PhotonTools.cc,v 1.12 2011/08/03 17:15:44 bendavid Exp $

#include "MitPhysics/Utils/interface/PhotonTools.h"
#include "MitPhysics/Utils/interface/ElectronTools.h"
#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitAna/DataTree/interface/StableData.h"
#include <TFile.h>
#include <TRandom3.h>

ClassImp(mithep::PhotonTools)

using namespace mithep;

//--------------------------------------------------------------------------------------------------
PhotonTools::PhotonTools()  
{
  // Constructor.
}


void PhotonTools::ScalePhoton(Photon* p, Double_t scale) {
  if( !p ) return;
  FourVectorM mom = p->Mom();
  p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());
  
}

void PhotonTools::SmearPhoton(Photon* p, TRandom3* rng, Double_t width, UInt_t iSeed) {
  
  if( !p  )       return;
  if( !rng)       return;
  if( width < 0.) return;

  rng->SetSeed(iSeed);
  FourVectorM mom = p->Mom();
  Double_t scale = rng->Gaus(1.,width);

  if( scale > 0)
    p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());
  
  return;
}

void PhotonTools::SmearPhotonError(Photon* p, Double_t width) {
  
  if( !p  )       return;
  if( width < 0.) return;  
  
  Double_t err = p->EnergyErrSmeared();
  if (err>=0.0) {
    p->SetEnergyErrSmeared(TMath::Sqrt(err*err + width*width*p->E()*p->E()));
  }
  
}

//--------------------------------------------------------------------------------------------------
Bool_t PhotonTools::PassConversionId(const Photon *p, const DecayParticle *c) {

  if (!c) return kTRUE;
  
  ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
  Double_t deta = c->Eta()-dirconvsc.Eta();
  Double_t dphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
  Double_t eoverp = p->SCluster()->Energy()/c->P();
  
  if (p->IsEB() && eoverp>2.0) return kFALSE;
  if (p->IsEE() && eoverp>3.0) return kFALSE;
  
  if (p->IsEE() && TMath::Abs(deta)>0.01) return kFALSE;
  if (p->IsEE() && TMath::Abs(dphi)>0.01) return kFALSE;

  return kTRUE;
    
}

//--------------------------------------------------------------------------------------------------
Bool_t PhotonTools::PassElectronVeto(const Photon *p, const ElectronCol *els) {

  Bool_t pass = kTRUE;
  for (UInt_t i=0; i<els->GetEntries(); ++i) {
    const Electron *e = els->At(i);
    if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0) {
      pass = kFALSE;
    }
  }
  
  return pass;
}

//--------------------------------------------------------------------------------------------------
const Electron *PhotonTools::MatchedElectron(const Photon *p, const ElectronCol *els) {

  for (UInt_t i=0; i<els->GetEntries(); ++i) {
    const Electron *e = els->At(i);
    if ( e->SCluster()==p->SCluster() ) {
      return e;
    }
  }
  
  return 0;
}

//--------------------------------------------------------------------------------------------------
const Photon *PhotonTools::MatchedPhoton(const Electron *e, const PhotonCol *phs) {

  for (UInt_t i=0; i<phs->GetEntries(); ++i) {
    const Photon *p = phs->At(i);
    if ( p->SCluster()==e->SCluster() ) {
      return p;
    }
  }
  
  return 0;
}

//--------------------------------------------------------------------------------------------------
const SuperCluster *PhotonTools::MatchedSC(const SuperCluster *psc, const SuperClusterCol *scs, Double_t drMin) {

  Double_t drsmallest = 999.;
  const SuperCluster *match = 0;
  for (UInt_t i=0; i<scs->GetEntries(); ++i) {
    const SuperCluster *sc = scs->At(i);
    Double_t dr = MathUtils::DeltaR(*sc,*psc);
    if ( dr<drsmallest && dr<drMin ) {
      drsmallest = dr;
      match = sc;
    }
  }
  
  return match;
}

//--------------------------------------------------------------------------------------------------
Double_t PhotonTools::ElectronVetoCiC(const Photon *p, const ElectronCol *els) {
  
  for (UInt_t i=0; i<els->GetEntries(); ++i) {
    const Electron *e = els->At(i);
    if (e->SCluster()==p->SCluster() ) {
      //if( e->GsfTrk()->NExpectedHitsInner()==0 && e->GsfTrk()->Pt() > 2.5 ) {
      if( e->GsfTrk()->NExpectedHitsInner()==0 ) {
        double dEta = e->DeltaEtaSuperClusterTrackAtVtx();
        double dPhi = e->DeltaPhiSuperClusterTrackAtVtx();
        double dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi);
        return dR;
      }
    }
  }  
  return 99.;
}

//--------------------------------------------------------------------------------------------------
Bool_t PhotonTools::PassElectronVetoConvRecovery(const Photon *p, const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *v) {

  Bool_t pass = kTRUE;
  for (UInt_t i=0; i<els->GetEntries(); ++i) {
    const Electron *e = els->At(i);
    if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0 && ElectronTools::PassConversionFilter(e, conversions, 
                                                         v, 0, 1e-6, 2.0, kFALSE, kFALSE) ) {
      pass = kFALSE;
    }
  }
  
  return pass;
}

//--------------------------------------------------------------------------------------------------
Bool_t PhotonTools::PassTriggerMatching(const Photon *p, const TriggerObjectCol *trigobjs)
{
  
  for (UInt_t i=0; i<trigobjs->GetEntries(); ++i) {
    const TriggerObject *trigobj = trigobjs->At(i);
    if (trigobj->TriggerType()==TriggerObject::TriggerCluster || trigobj->TriggerType()==TriggerObject::TriggerElectron || trigobj->TriggerType()==TriggerObject::TriggerPhoton) {
      if (MathUtils::DeltaR(p->SCluster(),trigobj)<0.3) {
        return kTRUE;
      }
    }
  }
  
  return kFALSE;
  
  
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedConversion(const Photon *p, const DecayParticleCol *conversions, 
                                               const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
                                               Double_t lxyMin, Double_t dRMin) {
  
  return MatchedConversion(p->SCluster(), conversions, vtx, nWrongHitsMax, probMin, lxyMin, dRMin);
  
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedConversion(const SuperCluster *sc, const DecayParticleCol *conversions, 
                                               const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
                                               Double_t lxyMin, Double_t dRMin) {
  
  const DecayParticle *match = 0;
  Double_t rhosmallest = 999.;
  for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
    const DecayParticle *c = conversions->At(i);
    ThreeVector dirconvsc = ThreeVector(sc->Point()) - c->Position();
    Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
    Double_t rho = c->Position().Rho();
    if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
      Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
      Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
      if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
        rhosmallest = rho;
        match = c;
      }
    }
    
  }
  
  return match;
  
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedConversion(const Track *t, const DecayParticleCol *conversions, 
                                               const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
                                               Double_t lxyMin) {
  
  for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
    const DecayParticle *c = conversions->At(i);
    if (c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
      Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
      Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
      if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
        const Track *ct1 = dynamic_cast<const ChargedParticle*>(c->Daughter(0))->Trk();
        const Track *ct2 = dynamic_cast<const ChargedParticle*>(c->Daughter(1))->Trk();
        if (t==ct1 || t==ct2) return c;
      }
    }
    
  }
  
  return 0;
  
}

PhotonTools::DiphotonR9EtaCats PhotonTools::DiphotonR9EtaCat(const Photon *p1, const Photon *p2) {
  
  if (p1->IsEB() && p2->IsEB()) {
    if (p1->R9()>0.93 && p2->R9()>0.93) return kCat1;
    else return kCat2;
    
  }
  else {
    if (p1->R9()>0.93 && p2->R9()>0.93) return kCat3;
    else return kCat4; 
  }
  
}

PhotonTools::DiphotonR9EtaConversionCats PhotonTools::DiphotonR9EtaConversionCat(const Photon *p1, const Photon *p2, const DecayParticleCol *conversions, const BaseVertex *v) {
  
  const DecayParticle *conv1 = MatchedConversion(p1, conversions, v);
  const DecayParticle *conv2 = MatchedConversion(p2, conversions, v);
    
  if (p1->IsEB() && p2->IsEB()) {
    if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat1;
    else if (conv1||conv2) return kNewCat2;
    else return kNewCat3;
    
  }
  else {
    if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat4;
    else if (conv1||conv2) return kNewCat5;
    else return kNewCat6;
  }
  
}

PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
  if( p->SCluster()->AbsEta()<1.5 ) {
    if ( p->R9() > 0.94 ) return kCiCCat1;
    else return kCiCCat2;
  } else {
    if ( p->R9() > 0.94 ) return kCiCCat3;
    else return kCiCCat4;
  }

  // shoud NEVER happen, but you never know...
  return kCiCNoCat;
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedCiCConversion(const Photon *p, const DecayParticleCol *conversions, 
                                                       Double_t dPhiMin,
                                                       Double_t dEtaMin,
                                                       Double_t dRMin,
                                                       bool print) {
  
  // if there are no conversons, return
  if ( !p || !conversions)  return NULL;

  const DecayParticle *match = NULL;

  double minDeta = 999.;
  double minDphi = 999.;
  double minDR   = 999.;

  double phPhi = p->SCluster()->Phi();
  double phEta = p->SCluster()->Eta();

  if(print)
    std::cout<<"  --- conv match photon eta = "<<phEta<<"  phi = "<<phPhi<<std::endl;  


  for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
    const DecayParticle *c = conversions->At(i);

    if(print)
      std::cout<< "   c "<<i+1<<"  pt = "<<c->Pt()<<std::endl;

    if(c->Pt()   < 1. )    continue; // is this refittedPirMomentum?

    //ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
    //ThreeVector dirconvsc = p->CaloPos() - c->Position();

    double convPhi = c->Phi();
    double convEta = c->Eta();
    
    const ThreeVector wrong(0,0,0);
    double Zvertex = c->DzCorrected(wrong);
    // ------------------ FROM GLOBE ----------------------
    //---Definitions for ECAL
    const float R_ECAL           = 136.5;
    const float Z_Endcap         = 328.0;
    const float etaBarrelEndcap  = 1.479; 
    
    //---ETA correction
    float Theta = 0.0  ; 
    float ZEcal = R_ECAL*sinh(convEta)+Zvertex;
    
    if(ZEcal != 0.0) Theta = TMath::ATan(R_ECAL/ZEcal);
    if(Theta<0.0) Theta = Theta+M_PI;
    double fEta = - TMath::Log(TMath::Tan(0.5*Theta));
    
    if( fabs(fEta) > etaBarrelEndcap ) {
      float Zend = Z_Endcap ;
      if(convEta<0.0 )  Zend = -Zend ;
      float Zlen = Zend - Zvertex ;
      float RR = Zlen/TMath::SinH(convEta); 
      Theta = TMath::ATan(RR/Zend);
      if(Theta<0.0) Theta = Theta+M_PI;
      fEta = -TMath::Log(TMath::Tan(0.5*Theta));              
    } 
    // ---------------------------------------------------

    if(print) {
      std::cout<<"          eta bare = "<<convEta<<std::endl;
      std::cout<<"          eta new  = "<<fEta<<std::endl;
      std::cout<<"          phi      = "<<convPhi<<std::endl;
    }

    convEta = fEta;
    
    Double_t dphi = TMath::Abs(phPhi - convPhi);
    if(dphi > M_PI) dphi = 2.*M_PI-dphi;
    //Double_t deta = c->Eta()-dirconvsc.Eta();
    Double_t deta = TMath::Abs(convEta-phEta);
    Double_t dR = TMath::Sqrt(dphi*dphi+deta*deta);
    //if(dphi < minDphi && TMath::Abs(deta)<minDeta) {
    if(dR < minDR) {
      minDR = dR;
      minDphi = dphi;
      minDeta = TMath::Abs(deta);
      match = c;

      if(print)
        std::cout<<" conv "<<i+1<<" matches with dPhi = "<<minDphi<<"   dEta = "<<minDeta<<std::endl;

    } 
  }
  
  //if(minDphi < dPhiMin && minDeta < dEtaMin)
  if(minDR < dRMin)
    return match;
  else 
    return NULL;
}

bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx, 
                                   const TrackCol* trackCol,
                                   const ElectronCol* eleCol,
                                   const VertexCol* vtxCol,
                                   double rho, double ptmin, 
                                   bool applyEleVeto,
                                   bool print, float* kin) {

  
  // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
  float cic4_allcuts_temp_sublead[] = { 
    3.8,         2.2,         1.77,        1.29,
    11.7,        3.4,         3.9,         1.84,
    3.5,         2.2,         2.3,         1.45,
    0.0106,      0.0097,      0.028,       0.027,
    0.082,       0.062,       0.065,       0.048,
    0.94,        0.36,        0.94,        0.32,
    1.,          0.062,       0.97,        0.97,
    1.5,         1.5,         1.5,         1.5 };  // the last line is PixelmatchVeto and un-used
  
  // cut on Et instead of Pt???    
  Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
    
  // compute all relevant observables first
  double ecalIso3 = ph->EcalRecHitIsoDr03();
  double ecalIso4 = ph->EcalRecHitIsoDr04();
  double hcalIso4 = ph->HcalTowerSumEtDr04();

  unsigned int wVtxInd = 0;

  double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);

  // track iso only
  double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);

  // track iso worst vtx
  double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
  
  double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
  double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
  
  double tIso1 = (combIso1) *50./ph->Et();
  double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
  //double tIso2 = (combIso2) *50./ph->Et();
  double tIso3 = (trackIso3)*50./ph->Et();
  
  double covIEtaIEta  =ph->CoviEtaiEta();
  double HoE = ph->HadOverEm();
  
  double R9 = ph->R9();
  
  double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);

  // check which category it is ...
  int _tCat = 1;
  if ( !isbarrel ) _tCat = 3;
  if ( R9 < 0.94 ) _tCat++;
  
  if(print) {
    std::cout<<" -------------------------- "<<std::endl;
    std::cout<<" photon Et  = "<<ph->Et()<<std::endl;
    std::cout<<"        Eta = "<<ph->SCluster()->Eta()<<std::endl;
    std::cout<<"        HoE = "<<HoE<<std::endl;
    std::cout<<"         R9 = "<<R9<<std::endl;
    std::cout<<"         dR = "<<dRTrack<<std::endl;
    std::cout<<"       iso1 = "<<tIso1<<std::endl; 
    std::cout<<"       iso2 = "<<tIso2<<std::endl; 
    std::cout<<"       iso3 = "<<tIso3<<std::endl; 
  }

  if(kin) {
    kin[0] = tIso1;
    kin[1] = tIso2;
    kin[2] = tIso3;
    kin[3] = covIEtaIEta;
    kin[4] = HoE;
    kin[5] = R9;
    kin[6] = dRTrack;
    kin[7] = (float) ph->Pt();
    kin[8] = (float) ph->Eta();
    kin[9] = (float) ph->Phi();

    // iso quantities separate
    kin[10] = ecalIso3;
    kin[11] = ecalIso4;
    kin[12] = hcalIso4;
    kin[13] = trackIso1;
    kin[14] = trackIso2;
    kin[15] = trackIso3;

    kin[16] = (float) ph->Et();
    kin[17] = (float) ph->E();

    kin[18] = (float) ph->SCluster()->Eta();
    kin[19] = (float) ph->SCluster()->Phi();
  }

  float passCuts = 1.;

  if ( ph->Pt()     <= ptmin      ) passCuts = -1.;

  // not needed anymore, do in pre-selection...
  //if (  ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
  
  if(   ! (    tIso1                          < cic4_allcuts_temp_sublead[_tCat-1+0*4]
               && tIso2                       < cic4_allcuts_temp_sublead[_tCat-1+1*4]
               && tIso3                       < cic4_allcuts_temp_sublead[_tCat-1+2*4]
               && covIEtaIEta                 < cic4_allcuts_temp_sublead[_tCat-1+3*4]
               && HoE                         < cic4_allcuts_temp_sublead[_tCat-1+4*4]
               && R9                          > cic4_allcuts_temp_sublead[_tCat-1+5*4]
               && ( dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] || !applyEleVeto ) ) )   passCuts = -1.;
  
  if(kin) {    
    kin[20] = passCuts;
    kin[21] = (float) _tCat;
  }    

  if(passCuts > 0.) return true;
  return false;
}

const MCParticle *PhotonTools::MatchMC(const Particle *ph, const MCParticleCol *c, Bool_t matchElectrons) {

//  printf("Start loop\n");
  for (UInt_t i=0; i<c->GetEntries(); ++i) {
    const MCParticle *p = c->At(i);
//     if (p->IsGenerated() && p->AbsPdgId()==11 && (p->DistinctMother()->AbsPdgId()==23|| p->DistinctMother()->AbsPdgId()==24)) {
//       printf("pdgid = %i, status = %i, pt = %5f, mother pdgid = %i\n",p->PdgId(),p->Status(),p->Pt(),p->Mother()->PdgId());
//     }
    if (matchElectrons && p->AbsPdgId()==11 && p->IsGenerated() && p->Status()==3 && MathUtils::DeltaR(*ph,*p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==23 || p->Mother()->AbsPdgId()==24 || p->Mother()->AbsPdgId()==22)) {
      return p;
    }
    if ( p->AbsPdgId()==22 && p->IsGenerated() && MathUtils::DeltaR(*ph,*p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==25 || p->Mother()->AbsPdgId()<=21) ) {
      return p;
    }
  }
  return 0;
}

PhotonTools::DiphotonR9EtaPtCats PhotonTools::DiphotonR9EtaPtCat(const Photon *p1, const Photon *p2) {

  PhotonTools::CiCBaseLineCats cat1 = CiCBaseLineCat(p1);
  PhotonTools::CiCBaseLineCats cat2 = CiCBaseLineCat(p2);
  
  PhotonTools::DiphotonR9EtaPtCats evtcat = PhotonTools::kOctCat0;
  
  bool ph1IsEB = (cat1 ==  PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
  bool ph2IsEB = (cat2 ==  PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);

  bool ph1IsHR9 = (cat1 ==  PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
  bool ph2IsHR9 = (cat2 ==  PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
  
  if( ph1IsEB && ph2IsEB )
    evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat0 : PhotonTools::kOctCat1);
  else
    evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat2 : PhotonTools::kOctCat3);
  
  float ptgg = (p1->Mom()+p2->Mom()).Pt();
  if (ptgg<40.0) evtcat = PhotonTools::DiphotonR9EtaPtCats(evtcat + 4);
  
  
  return evtcat;
}
Revision:	1.13
Committed:	Thu Sep 8 15:51:24 2011 UTC (13 years, 7 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025pre2, Mit_025pre1
Changes since 1.12:	+72 -7 lines
Log Message:	Add tool and required changes for photon energy regression
#	User	Rev	Content
1	bendavid	1.13	// $Id: PhotonTools.cc,v 1.12 2011/08/03 17:15:44 bendavid Exp $
2	bendavid	1.1
3			#include "MitPhysics/Utils/interface/PhotonTools.h"
4			#include "MitPhysics/Utils/interface/ElectronTools.h"
5	fabstoec	1.5	#include "MitPhysics/Utils/interface/IsolationTools.h"
6	bendavid	1.1	#include "MitAna/DataTree/interface/StableData.h"
7			#include <TFile.h>
8	fabstoec	1.8	#include <TRandom3.h>
9	bendavid	1.1
10			ClassImp(mithep::PhotonTools)
11
12			using namespace mithep;
13
14			//--------------------------------------------------------------------------------------------------
15			PhotonTools::PhotonTools()
16			{
17			// Constructor.
18			}
19
20	fabstoec	1.8
21			void PhotonTools::ScalePhoton(Photon* p, Double_t scale) {
22			if( !p ) return;
23			FourVectorM mom = p->Mom();
24			p->SetMom(scalemom.X(), scalemom.Y(), scalemom.Z(), scalemom.E());
25
26			}
27
28			void PhotonTools::SmearPhoton(Photon* p, TRandom3* rng, Double_t width, UInt_t iSeed) {
29
30			if( !p ) return;
31			if( !rng) return;
32			if( width < 0.) return;
33
34	fabstoec	1.10	rng->SetSeed(iSeed);
35	fabstoec	1.8	FourVectorM mom = p->Mom();
36			Double_t scale = rng->Gaus(1.,width);
37	fabstoec	1.10
38	fabstoec	1.8	if( scale > 0)
39			p->SetMom(scalemom.X(), scalemom.Y(), scalemom.Z(), scalemom.E());
40	bendavid	1.13
41			return;
42			}
43	fabstoec	1.8
44	bendavid	1.13	void PhotonTools::SmearPhotonError(Photon* p, Double_t width) {
45
46			if( !p ) return;
47			if( width < 0.) return;
48
49			Double_t err = p->EnergyErrSmeared();
50			if (err>=0.0) {
51			p->SetEnergyErrSmeared(TMath::Sqrt(errerr + widthwidthp->E()p->E()));
52			}
53
54	fabstoec	1.8	}
55
56	bendavid	1.1	//--------------------------------------------------------------------------------------------------
57			Bool_t PhotonTools::PassConversionId(const Photon p, const DecayParticle c) {
58
59			if (!c) return kTRUE;
60
61			ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
62			Double_t deta = c->Eta()-dirconvsc.Eta();
63			Double_t dphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
64			Double_t eoverp = p->SCluster()->Energy()/c->P();
65
66			if (p->IsEB() && eoverp>2.0) return kFALSE;
67			if (p->IsEE() && eoverp>3.0) return kFALSE;
68
69			if (p->IsEE() && TMath::Abs(deta)>0.01) return kFALSE;
70			if (p->IsEE() && TMath::Abs(dphi)>0.01) return kFALSE;
71
72			return kTRUE;
73
74			}
75
76			//--------------------------------------------------------------------------------------------------
77			Bool_t PhotonTools::PassElectronVeto(const Photon p, const ElectronCol els) {
78
79			Bool_t pass = kTRUE;
80			for (UInt_t i=0; i<els->GetEntries(); ++i) {
81			const Electron *e = els->At(i);
82			if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0) {
83			pass = kFALSE;
84			}
85			}
86
87			return pass;
88			}
89
90			//--------------------------------------------------------------------------------------------------
91	bendavid	1.13	const Electron PhotonTools::MatchedElectron(const Photon p, const ElectronCol *els) {
92
93			for (UInt_t i=0; i<els->GetEntries(); ++i) {
94			const Electron *e = els->At(i);
95			if ( e->SCluster()==p->SCluster() ) {
96			return e;
97			}
98			}
99
100			return 0;
101			}
102
103			//--------------------------------------------------------------------------------------------------
104			const Photon PhotonTools::MatchedPhoton(const Electron e, const PhotonCol *phs) {
105
106			for (UInt_t i=0; i<phs->GetEntries(); ++i) {
107			const Photon *p = phs->At(i);
108			if ( p->SCluster()==e->SCluster() ) {
109			return p;
110			}
111			}
112
113			return 0;
114			}
115
116			//--------------------------------------------------------------------------------------------------
117			const SuperCluster PhotonTools::MatchedSC(const SuperCluster psc, const SuperClusterCol *scs, Double_t drMin) {
118
119			Double_t drsmallest = 999.;
120			const SuperCluster *match = 0;
121			for (UInt_t i=0; i<scs->GetEntries(); ++i) {
122			const SuperCluster *sc = scs->At(i);
123			Double_t dr = MathUtils::DeltaR(sc,psc);
124			if ( dr<drsmallest && dr<drMin ) {
125			drsmallest = dr;
126			match = sc;
127			}
128			}
129
130			return match;
131			}
132
133			//--------------------------------------------------------------------------------------------------
134	fabstoec	1.4	Double_t PhotonTools::ElectronVetoCiC(const Photon p, const ElectronCol els) {
135
136			for (UInt_t i=0; i<els->GetEntries(); ++i) {
137			const Electron *e = els->At(i);
138	fabstoec	1.5	if (e->SCluster()==p->SCluster() ) {
139			//if( e->GsfTrk()->NExpectedHitsInner()==0 && e->GsfTrk()->Pt() > 2.5 ) {
140			if( e->GsfTrk()->NExpectedHitsInner()==0 ) {
141			double dEta = e->DeltaEtaSuperClusterTrackAtVtx();
142			double dPhi = e->DeltaPhiSuperClusterTrackAtVtx();
143			double dR = TMath::Sqrt(dEtadEta+dPhidPhi);
144			return dR;
145			}
146			}
147			}
148			return 99.;
149	fabstoec	1.4	}
150
151			//--------------------------------------------------------------------------------------------------
152	bendavid	1.1	Bool_t PhotonTools::PassElectronVetoConvRecovery(const Photon p, const ElectronCol els, const DecayParticleCol conversions, const BaseVertex v) {
153
154			Bool_t pass = kTRUE;
155			for (UInt_t i=0; i<els->GetEntries(); ++i) {
156			const Electron *e = els->At(i);
157			if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0 && ElectronTools::PassConversionFilter(e, conversions,
158			v, 0, 1e-6, 2.0, kFALSE, kFALSE) ) {
159			pass = kFALSE;
160			}
161			}
162
163			return pass;
164			}
165
166			//--------------------------------------------------------------------------------------------------
167			Bool_t PhotonTools::PassTriggerMatching(const Photon p, const TriggerObjectCol trigobjs)
168			{
169
170			for (UInt_t i=0; i<trigobjs->GetEntries(); ++i) {
171			const TriggerObject *trigobj = trigobjs->At(i);
172			if (trigobj->TriggerType()==TriggerObject::TriggerCluster \|\| trigobj->TriggerType()==TriggerObject::TriggerElectron \|\| trigobj->TriggerType()==TriggerObject::TriggerPhoton) {
173			if (MathUtils::DeltaR(p->SCluster(),trigobj)<0.3) {
174			return kTRUE;
175			}
176			}
177			}
178
179			return kFALSE;
180
181
182			}
183
184			//--------------------------------------------------------------------------------------------------
185			const DecayParticle PhotonTools::MatchedConversion(const Photon p, const DecayParticleCol *conversions,
186			const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
187			Double_t lxyMin, Double_t dRMin) {
188
189	bendavid	1.13	return MatchedConversion(p->SCluster(), conversions, vtx, nWrongHitsMax, probMin, lxyMin, dRMin);
190
191			}
192
193			//--------------------------------------------------------------------------------------------------
194			const DecayParticle PhotonTools::MatchedConversion(const SuperCluster sc, const DecayParticleCol *conversions,
195			const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
196			Double_t lxyMin, Double_t dRMin) {
197
198	bendavid	1.1	const DecayParticle *match = 0;
199	bendavid	1.3	Double_t rhosmallest = 999.;
200	bendavid	1.1	for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
201			const DecayParticle *c = conversions->At(i);
202	bendavid	1.13	ThreeVector dirconvsc = ThreeVector(sc->Point()) - c->Position();
203	bendavid	1.1	Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
204	bendavid	1.3	Double_t rho = c->Position().Rho();
205			if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
206	bendavid	1.1	Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
207			Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
208			if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
209	bendavid	1.3	rhosmallest = rho;
210	bendavid	1.1	match = c;
211			}
212			}
213
214			}
215
216			return match;
217
218			}
219
220	bendavid	1.2	//--------------------------------------------------------------------------------------------------
221			const DecayParticle PhotonTools::MatchedConversion(const Track t, const DecayParticleCol *conversions,
222			const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
223			Double_t lxyMin) {
224
225			for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
226			const DecayParticle *c = conversions->At(i);
227			if (c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
228			Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
229			Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
230			if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
231			const Track ct1 = dynamic_cast<const ChargedParticle>(c->Daughter(0))->Trk();
232			const Track ct2 = dynamic_cast<const ChargedParticle>(c->Daughter(1))->Trk();
233			if (t==ct1 \|\| t==ct2) return c;
234			}
235			}
236
237			}
238
239			return 0;
240
241			}
242
243	bendavid	1.1	PhotonTools::DiphotonR9EtaCats PhotonTools::DiphotonR9EtaCat(const Photon p1, const Photon p2) {
244
245			if (p1->IsEB() && p2->IsEB()) {
246			if (p1->R9()>0.93 && p2->R9()>0.93) return kCat1;
247			else return kCat2;
248
249			}
250			else {
251			if (p1->R9()>0.93 && p2->R9()>0.93) return kCat3;
252			else return kCat4;
253			}
254
255			}
256
257			PhotonTools::DiphotonR9EtaConversionCats PhotonTools::DiphotonR9EtaConversionCat(const Photon p1, const Photon p2, const DecayParticleCol conversions, const BaseVertex v) {
258
259			const DecayParticle *conv1 = MatchedConversion(p1, conversions, v);
260			const DecayParticle *conv2 = MatchedConversion(p2, conversions, v);
261
262			if (p1->IsEB() && p2->IsEB()) {
263			if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat1;
264			else if (conv1\|\|conv2) return kNewCat2;
265			else return kNewCat3;
266
267			}
268			else {
269			if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat4;
270			else if (conv1\|\|conv2) return kNewCat5;
271			else return kNewCat6;
272			}
273
274	fabstoec	1.4	}
275
276			PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
277	bendavid	1.12	if( p->SCluster()->AbsEta()<1.5 ) {
278	fabstoec	1.4	if ( p->R9() > 0.94 ) return kCiCCat1;
279			else return kCiCCat2;
280			} else {
281			if ( p->R9() > 0.94 ) return kCiCCat3;
282			else return kCiCCat4;
283			}
284	fabstoec	1.8
285			// shoud NEVER happen, but you never know...
286			return kCiCNoCat;
287	fabstoec	1.4	}
288	fabstoec	1.5
289			//--------------------------------------------------------------------------------------------------
290			const DecayParticle PhotonTools::MatchedCiCConversion(const Photon p, const DecayParticleCol *conversions,
291			Double_t dPhiMin,
292	fabstoec	1.7	Double_t dEtaMin,
293			Double_t dRMin,
294			bool print) {
295	fabstoec	1.5
296	fabstoec	1.8	// if there are no conversons, return
297			if ( !p \|\| !conversions) return NULL;
298
299	fabstoec	1.6	const DecayParticle *match = NULL;
300
301			double minDeta = 999.;
302			double minDphi = 999.;
303	fabstoec	1.7	double minDR = 999.;
304	fabstoec	1.6
305			double phPhi = p->SCluster()->Phi();
306			double phEta = p->SCluster()->Eta();
307	fabstoec	1.7
308			if(print)
309			std::cout<<" --- conv match photon eta = "<<phEta<<" phi = "<<phPhi<<std::endl;
310
311
312	fabstoec	1.5	for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
313			const DecayParticle *c = conversions->At(i);
314	fabstoec	1.6
315	fabstoec	1.7	if(print)
316			std::cout<< " c "<<i+1<<" pt = "<<c->Pt()<<std::endl;
317
318	fabstoec	1.6	if(c->Pt() < 1. ) continue; // is this refittedPirMomentum?
319
320			//ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
321			//ThreeVector dirconvsc = p->CaloPos() - c->Position();
322
323			double convPhi = c->Phi();
324			double convEta = c->Eta();
325	fabstoec	1.7
326	fabstoec	1.6	const ThreeVector wrong(0,0,0);
327			double Zvertex = c->DzCorrected(wrong);
328			// ------------------ FROM GLOBE ----------------------
329			//---Definitions for ECAL
330			const float R_ECAL = 136.5;
331			const float Z_Endcap = 328.0;
332			const float etaBarrelEndcap = 1.479;
333	fabstoec	1.7
334	fabstoec	1.6	//---ETA correction
335			float Theta = 0.0 ;
336			float ZEcal = R_ECAL*sinh(convEta)+Zvertex;
337	fabstoec	1.7
338	fabstoec	1.6	if(ZEcal != 0.0) Theta = TMath::ATan(R_ECAL/ZEcal);
339			if(Theta<0.0) Theta = Theta+M_PI;
340			double fEta = - TMath::Log(TMath::Tan(0.5*Theta));
341	fabstoec	1.5
342	fabstoec	1.6	if( fabs(fEta) > etaBarrelEndcap ) {
343			float Zend = Z_Endcap ;
344			if(convEta<0.0 ) Zend = -Zend ;
345			float Zlen = Zend - Zvertex ;
346			float RR = Zlen/TMath::SinH(convEta);
347			Theta = TMath::ATan(RR/Zend);
348			if(Theta<0.0) Theta = Theta+M_PI;
349			fEta = -TMath::Log(TMath::Tan(0.5*Theta));
350			}
351			// ---------------------------------------------------
352
353	fabstoec	1.7	if(print) {
354			std::cout<<" eta bare = "<<convEta<<std::endl;
355			std::cout<<" eta new = "<<fEta<<std::endl;
356			std::cout<<" phi = "<<convPhi<<std::endl;
357			}
358
359			convEta = fEta;
360
361	fabstoec	1.6	Double_t dphi = TMath::Abs(phPhi - convPhi);
362			if(dphi > M_PI) dphi = 2.*M_PI-dphi;
363			//Double_t deta = c->Eta()-dirconvsc.Eta();
364			Double_t deta = TMath::Abs(convEta-phEta);
365	fabstoec	1.7	Double_t dR = TMath::Sqrt(dphidphi+detadeta);
366			//if(dphi < minDphi && TMath::Abs(deta)<minDeta) {
367			if(dR < minDR) {
368			minDR = dR;
369			minDphi = dphi;
370			minDeta = TMath::Abs(deta);
371			match = c;
372
373			if(print)
374			std::cout<<" conv "<<i+1<<" matches with dPhi = "<<minDphi<<" dEta = "<<minDeta<<std::endl;
375
376	fabstoec	1.6	}
377	fabstoec	1.5	}
378
379	fabstoec	1.7	//if(minDphi < dPhiMin && minDeta < dEtaMin)
380			if(minDR < dRMin)
381			return match;
382			else
383			return NULL;
384	fabstoec	1.5	}
385
386	fabstoec	1.9	bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx,
387			const TrackCol* trackCol,
388			const ElectronCol* eleCol,
389			const VertexCol* vtxCol,
390	fabstoec	1.5	double rho, double ptmin,
391	fabstoec	1.9	bool applyEleVeto,
392	fabstoec	1.5	bool print, float* kin) {
393	fabstoec	1.9
394	fabstoec	1.5
395			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
396			float cic4_allcuts_temp_sublead[] = {
397			3.8, 2.2, 1.77, 1.29,
398			11.7, 3.4, 3.9, 1.84,
399			3.5, 2.2, 2.3, 1.45,
400			0.0106, 0.0097, 0.028, 0.027,
401			0.082, 0.062, 0.065, 0.048,
402			0.94, 0.36, 0.94, 0.32,
403			1., 0.062, 0.97, 0.97,
404			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
405
406			// cut on Et instead of Pt???
407			Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
408
409			// compute all relevant observables first
410			double ecalIso3 = ph->EcalRecHitIsoDr03();
411			double ecalIso4 = ph->EcalRecHitIsoDr04();
412			double hcalIso4 = ph->HcalTowerSumEtDr04();
413
414			unsigned int wVtxInd = 0;
415
416			double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
417
418			// track iso only
419			double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
420
421			// track iso worst vtx
422			double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
423
424			double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
425			double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
426
427			double tIso1 = (combIso1) *50./ph->Et();
428			double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
429			//double tIso2 = (combIso2) *50./ph->Et();
430			double tIso3 = (trackIso3)*50./ph->Et();
431
432			double covIEtaIEta =ph->CoviEtaiEta();
433			double HoE = ph->HadOverEm();
434
435			double R9 = ph->R9();
436
437			double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);
438
439			// check which category it is ...
440			int _tCat = 1;
441			if ( !isbarrel ) _tCat = 3;
442			if ( R9 < 0.94 ) _tCat++;
443
444			if(print) {
445			std::cout<<" -------------------------- "<<std::endl;
446			std::cout<<" photon Et = "<<ph->Et()<<std::endl;
447			std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
448			std::cout<<" HoE = "<<HoE<<std::endl;
449			std::cout<<" R9 = "<<R9<<std::endl;
450			std::cout<<" dR = "<<dRTrack<<std::endl;
451			std::cout<<" iso1 = "<<tIso1<<std::endl;
452			std::cout<<" iso2 = "<<tIso2<<std::endl;
453			std::cout<<" iso3 = "<<tIso3<<std::endl;
454			}
455
456			if(kin) {
457			kin[0] = tIso1;
458			kin[1] = tIso2;
459			kin[2] = tIso3;
460			kin[3] = covIEtaIEta;
461			kin[4] = HoE;
462			kin[5] = R9;
463			kin[6] = dRTrack;
464			kin[7] = (float) ph->Pt();
465			kin[8] = (float) ph->Eta();
466			kin[9] = (float) ph->Phi();
467
468			// iso quantities separate
469			kin[10] = ecalIso3;
470			kin[11] = ecalIso4;
471			kin[12] = hcalIso4;
472			kin[13] = trackIso1;
473			kin[14] = trackIso2;
474			kin[15] = trackIso3;
475
476			kin[16] = (float) ph->Et();
477			kin[17] = (float) ph->E();
478
479	bendavid	1.11	kin[18] = (float) ph->SCluster()->Eta();
480			kin[19] = (float) ph->SCluster()->Phi();
481	fabstoec	1.5	}
482
483			float passCuts = 1.;
484
485			if ( ph->Pt() <= ptmin ) passCuts = -1.;
486
487			// not needed anymore, do in pre-selection...
488			//if ( ph->SCluster()->AbsEta()>=2.5 \|\| (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
489
490			if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
491			&& tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
492			&& tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
493			&& covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
494			&& HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
495			&& R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4]
496	fabstoec	1.9	&& ( dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] \|\| !applyEleVeto ) ) ) passCuts = -1.;
497	fabstoec	1.5
498			if(kin) {
499	bendavid	1.11	kin[20] = passCuts;
500			kin[21] = (float) _tCat;
501	fabstoec	1.5	}
502
503			if(passCuts > 0.) return true;
504			return false;
505			}
506	bendavid	1.11
507	bendavid	1.13	const MCParticle PhotonTools::MatchMC(const Particle ph, const MCParticleCol *c, Bool_t matchElectrons) {
508
509			// printf("Start loop\n");
510	bendavid	1.11	for (UInt_t i=0; i<c->GetEntries(); ++i) {
511			const MCParticle *p = c->At(i);
512	bendavid	1.13	// if (p->IsGenerated() && p->AbsPdgId()==11 && (p->DistinctMother()->AbsPdgId()==23\|\| p->DistinctMother()->AbsPdgId()==24)) {
513			// printf("pdgid = %i, status = %i, pt = %5f, mother pdgid = %i\n",p->PdgId(),p->Status(),p->Pt(),p->Mother()->PdgId());
514	bendavid	1.12	// }
515			if (matchElectrons && p->AbsPdgId()==11 && p->IsGenerated() && p->Status()==3 && MathUtils::DeltaR(ph,p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==23 \|\| p->Mother()->AbsPdgId()==24 \|\| p->Mother()->AbsPdgId()==22)) {
516			return p;
517			}
518	bendavid	1.11	if ( p->AbsPdgId()==22 && p->IsGenerated() && MathUtils::DeltaR(ph,p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==25 \|\| p->Mother()->AbsPdgId()<=21) ) {
519			return p;
520			}
521			}
522			return 0;
523	bendavid	1.12	}
524
525			PhotonTools::DiphotonR9EtaPtCats PhotonTools::DiphotonR9EtaPtCat(const Photon p1, const Photon p2) {
526
527			PhotonTools::CiCBaseLineCats cat1 = CiCBaseLineCat(p1);
528			PhotonTools::CiCBaseLineCats cat2 = CiCBaseLineCat(p2);
529
530			PhotonTools::DiphotonR9EtaPtCats evtcat = PhotonTools::kOctCat0;
531
532			bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 \|\| cat1 == PhotonTools::kCiCCat2);
533			bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 \|\| cat2 == PhotonTools::kCiCCat2);
534
535			bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 \|\| cat1 == PhotonTools::kCiCCat3);
536			bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 \|\| cat2 == PhotonTools::kCiCCat3);
537
538			if( ph1IsEB && ph2IsEB )
539			evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat0 : PhotonTools::kOctCat1);
540			else
541			evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat2 : PhotonTools::kOctCat3);
542
543			float ptgg = (p1->Mom()+p2->Mom()).Pt();
544			if (ptgg<40.0) evtcat = PhotonTools::DiphotonR9EtaPtCats(evtcat + 4);
545
546
547			return evtcat;
548	bendavid	1.11	}