Mods/src/PhotonCiCMod.cc

#include "MitPhysics/Mods/interface/PhotonCiCMod.h"
#include "MitAna/DataTree/interface/PhotonCol.h"
#include "MitPhysics/Init/interface/ModNames.h"
#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitPhysics/Utils/interface/PhotonTools.h"
#include <TNtuple.h>

using namespace mithep;

ClassImp(mithep::PhotonCiCMod)

//--------------------------------------------------------------------------------------------------
PhotonCiCMod::PhotonCiCMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fPhotonBranchName  (Names::gkPhotonBrn),
  fGoodPhotonsName   (ModNames::gkGoodPhotonsName),
  fTrackBranchName   (Names::gkTrackBrn),
  fPileUpDenName     (Names::gkPileupEnergyDensityBrn),
  fElectronName      ("Electrons"),
  fPhotonPtMin(20.0),
  fApplySpikeRemoval(kFALSE),
  fAbsEtaMax(999.99),
  fPhotons(0),
  fTracks(0),
  fPileUpDen(0),
  fElectrons(0),
  
  fPVName(Names::gkPVBeamSpotBrn),
  fPV(0),
  fPVFromBranch(kTRUE),

  fConversions(0),
  fConversionName(Names::gkMvfConversionBrn),

  fBeamspot(0)

{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
void PhotonCiCMod::Process()
{
  // Process entries of the tree. 

  LoadEventObject(fPhotonBranchName,   fPhotons);

  PhotonOArr *GoodPhotons = new PhotonOArr;
  GoodPhotons->SetName(fGoodPhotonsName);
  GoodPhotons->SetOwner(kTRUE);

  Double_t _tRho = -1.;
  LoadEventObject(fTrackBranchName,    fTracks);
  LoadEventObject(fPileUpDenName,      fPileUpDen);
  LoadEventObject(fElectronName,       fElectrons);
  LoadEventObject(fPVName,             fPV);    
  LoadEventObject(fConversionName,     fConversions);
  LoadBranch("BeamSpot");
  
  if(fPileUpDen->GetEntries() > 0)
    _tRho = (Double_t) fPileUpDen->At(0)->Rho();    
  
  _tRho = 0.;

  bool doVtxSelection = true;

  const EventHeader* evtHead = this->GetEventHeader();

  unsigned int evtNum = evtHead->EvtNum();
  Float_t _evtNum1 = (Float_t) ( (int) (evtNum/10000.) );
  Float_t _evtNum2 = (Float_t) ( (int) (evtNum % 10000)  );

  double evtNumTest = (int) ( ( (double) _evtNum1 )*10000. + (double) _evtNum2 );

  Float_t _runNum  = (Float_t) evtHead->RunNum();
  Float_t _lumiSec = (Float_t) evtHead->LumiSec();


  unsigned int numVertices = fPV->GetEntries();

  const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
  
  PhotonOArr* preselPh  = new PhotonOArr;
  
  // 1. we do the pre-selection; but keep the non-passing photons in a secont container...
  for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {    
    const Photon *ph = fPhotons->At(i);

    if(ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) continue;
    if(ph->Et() < 20.)     continue;
    if(ph->HadOverEm() > 0.15) continue;
    if(ph->AbsEta() < 1.5) {
      if(ph->CoviEtaiEta() > 0.013) continue;
      //if(ph->EcalRecHitIsoDr03() > 10.) continue;
    } else {
      if(ph->CoviEtaiEta() > 0.03) continue;
      //if(ph->EcalRecHitIsoDr03() > 10.) continue;
    }

    Bool_t passSpikeRemovalFilter = kTRUE;
  
    if (ph->SCluster() && ph->SCluster()->Seed()) {
      if(ph->SCluster()->Seed()->Energy() > 5.0 && 
         ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
        passSpikeRemovalFilter = kFALSE; 
   }    
    if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;

    preselPh->Add(ph);
  }

  // sort both by pt... again ;)
  preselPh->Sort();

  unsigned int numPairs = 0;
  if( preselPh->GetEntries() > 0)
    numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
  
  // create all possible pairs of pre-selected photons

  std::vector<unsigned int> idxFst;
  std::vector<unsigned int> idxSec;
  std::vector<bool> pairPasses;

  if(numPairs > 0) {
    for(unsigned int iFst = 0; iFst <preselPh->GetEntries() - 1; ++iFst) {
      for(unsigned int iSec = iFst + 1; iSec <preselPh->GetEntries(); ++iSec) {
        idxFst.push_back(iFst);
        idxSec.push_back(iSec);
        pairPasses.push_back(true);
      }
    }
  }
  
  // array to store the index of 'chosen Vtx' for each pair
  const Vertex** theVtx = new const Vertex*[numPairs];
  
  // arays to store the Vtx 'fixed' photons
  Photon** fixPhFst = new Photon*[numPairs];
  Photon** fixPhSec = new Photon*[numPairs];

  // sotr pair-indices for pairs passing the selection
  std::vector<unsigned int> passPairs;
  passPairs.resize(0);

  unsigned int theChosenVtx = 0;

  float kinPh1[20];
  float kinPh2[20];

  for(int i=0; i<10;  ++i) {
    kinPh1[i] =-99.;
    kinPh2[i] =-99.;
  }

  float ptBefore1 = -99.;
  float ptBefore2 = -99.;

  bool print = false;
  //if(evtNum == 1677 || evtNum == 1943 || evtNum == 1694) {
  if(evtNum == 1943 && false) {
    std::cout<<" ------------------------------------------- "<<std::endl;
    std::cout<<"   printing info for event #"<<evtNum<<std::endl;
    print = true;
  }

  for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {    
    
    // copy the photons for manipulation
    fixPhFst[iPair] = new Photon(*preselPh->At(idxFst[iPair]));
    fixPhSec[iPair] = new Photon(*preselPh->At(idxSec[iPair]));

    // store the vertex for this pair (TODO: conversion Vtx)
    if(doVtxSelection) {
      unsigned int iVtx = findBestVertex(fixPhFst[iPair],fixPhSec[iPair],bsp, print);
      theVtx[iPair] =  fPV->At(iVtx);
      if(iPair == 0) theChosenVtx = iVtx;
    } else
      theVtx[iPair] =  fPV->At(0);

    if(iPair==0) {
      ptBefore1 = fixPhFst[iPair]->Pt();
      ptBefore2 = fixPhSec[iPair]->Pt();
    }

    // fix the kinematics for both events
    FourVectorM newMomFst = fixPhFst[iPair]->MomVtx(theVtx[iPair]->Position());
    FourVectorM newMomSec = fixPhSec[iPair]->MomVtx(theVtx[iPair]->Position());
    fixPhFst[iPair]->SetMom(newMomFst.X(), newMomFst.Y(), newMomFst.Z(), newMomFst.E());
    fixPhSec[iPair]->SetMom(newMomSec.X(), newMomSec.Y(), newMomSec.Z(), newMomSec.E());
    
    if(iPair != 0) {
      // check if both photons pass the CiC selection
      bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., false);
      bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., false);
      if( pass1 && pass2 )
        passPairs.push_back(iPair);
    } else {
      bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., false, kinPh1);
      bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., false, kinPh2);
      if( pass1 && pass2 )
        passPairs.push_back(iPair);
    }
  }
   

  // loop over all passing pairs and find the one with the highest sum Et
  Photon* phHard = NULL;
  Photon* phSoft = NULL;
  double maxSumEt = 0.;
  for(unsigned int iPair=0; iPair<passPairs.size(); ++iPair){
    double sumEt = fixPhFst[passPairs[iPair]]->Et();
    sumEt += fixPhSec[passPairs[iPair]]->Et();
    if( sumEt > maxSumEt ) {
      maxSumEt = sumEt;
      phHard = fixPhFst[passPairs[iPair]];
      phSoft = fixPhSec[passPairs[iPair]];
    }
  }

  if(phHard && phSoft) {
    GoodPhotons->AddOwned(phHard);
    GoodPhotons->AddOwned(phSoft);
  }

  // sort according to pt
  GoodPhotons->Sort();
  
  // add to event for other modules to use
  AddObjThisEvt(GoodPhotons);

  delete preselPh;


  bool doFill = (phHard && phSoft);
  Float_t _mass = ( doFill ? (phHard->Mom()+phSoft->Mom()).M() : -100.);
  Float_t _ptgg = ( doFill ? (phHard->Mom()+phSoft->Mom()).Pt() : -100.);


  Float_t fillEvent[] = { _mass,
                          _ptgg,
                          _evtNum1,
                          _evtNum2,
                          _runNum,
                          _lumiSec,
                          (float) theChosenVtx,
                          (float) numPairs,
                          kinPh1[0],
                          kinPh1[1],
                          kinPh1[2],
                          kinPh1[3],
                          kinPh1[4],
                          kinPh1[5],
                          kinPh1[6],
                          kinPh1[7],
                          kinPh1[8],
                          kinPh1[9],
                          kinPh1[10],
                          kinPh1[11],
                          kinPh1[12],
                          kinPh1[13],
                          kinPh1[14],
                          kinPh1[15],
                          kinPh1[16],
                          kinPh1[17],
                          kinPh1[18],
                          kinPh1[19],
                          kinPh2[0],
                          kinPh2[1],
                          kinPh2[2],
                          kinPh2[3],
                          kinPh2[4],
                          kinPh2[5],
                          kinPh2[6],
                          kinPh2[7],
                          kinPh2[8],
                          kinPh2[9],
                          kinPh2[10],
                          kinPh2[11],
                          kinPh2[12],
                          kinPh2[13],
                          kinPh2[14],
                          kinPh2[15],
                          kinPh2[16],
                          kinPh2[17],
                          kinPh2[18],
                          kinPh2[19],
                          ptBefore1,
                          ptBefore2
  };

  hCiCTuple->Fill(fillEvent);

  return;

}

//--------------------------------------------------------------------------------------------------
void PhotonCiCMod::SlaveBegin()
{
  // Run startup code on the computer (slave) doing the actual analysis. Here,
  // we just request the photon collection branch.

  ReqEventObject(fPhotonBranchName,   fPhotons,   kTRUE);
  ReqEventObject(fTrackBranchName,    fTracks,    kTRUE);
  ReqEventObject(fElectronName,       fElectrons, kTRUE);  
  ReqEventObject(fPileUpDenName,      fPileUpDen, kTRUE);
  ReqEventObject(fPVName,             fPV,        fPVFromBranch);
  ReqEventObject(fConversionName,     fConversions,kTRUE);
  ReqBranch("BeamSpot",fBeamspot);


  hCiCTuple = new TNtuple("hCiCTuple","hCiCTuple","mass:ptgg:evtnum1:evtnum2:runnum:lumisec:ivtx:npairs:ph1Iso1:ph1Iso2:ph1Iso3:ph1Cov:ph1HoE:ph1R9:ph1DR:ph1Pt:ph1Eta:ph1Phi:ph1Eiso3:ph1Eiso4:ph1Hiso4:ph1TisoA:ph1TisoW:ph1Tiso:ph1Et:ph1E:ph1Pass:ph1Cat:ph2Iso1:ph2Iso2:ph2Iso3:ph2Cov:ph2HoE:ph2R9:ph2DR:ph2Pt:ph2Eta:ph2Phi:ph2Eiso3:ph2Eiso4:ph2Hiso4:ph2TisoA:ph2TisoW:ph2Tiso:ph2Et:ph2E:ph2Pass:ph2Cat:ph1UPt:ph2UPt");
  
  AddOutput(hCiCTuple);

}

// return the index of the bext vertex
unsigned int PhotonCiCMod::findBestVertex(Photon* ph1, Photon* ph2, const BaseVertex *bsp, bool print) {
  
  // loop over all vertices and assigne the ranks
  int* ptbal_rank  = new int[fPV->GetEntries()];
  int* ptasym_rank = new int[fPV->GetEntries()];
  int* total_rank  = new int[fPV->GetEntries()];
  double* ptbal = new double[fPV->GetEntries()];
  double* ptasym = new double[fPV->GetEntries()];
  
  unsigned int numVertices = fPV->GetEntries();

  double ptgg = 0.;    // stored for later in the conversion

  if(print) {
    std::cout<<" --------------------------------- "<<std::endl;
    std::cout<<"   looking for Vtx for photon pair "<<std::endl;
    std::cout<<"                            pt 1 = "<<ph1->Et()<<std::endl;
    std::cout<<"                            pt 2 = "<<ph2->Et()<<std::endl;
    std::cout<<"                    among #"<<numVertices<<" Vertices."<<std::endl;
  }


  for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
    if(print)
      std::cout<<std::endl<<"       Vertex #"<<iVtx<<std::endl;

    const Vertex* tVtx = fPV->At(iVtx);
    ptbal [iVtx] = 0.0;
    ptasym[iVtx] = 0.0;
    ptbal_rank [iVtx] = 1;
    ptasym_rank[iVtx] = 1;
    
    // compute the photon momenta with respect to this Vtx
    FourVectorM newMomFst = ph1->MomVtx(tVtx->Position());
    FourVectorM newMomSec = ph2->MomVtx(tVtx->Position());
    
    FourVectorM higgsMom = newMomFst+newMomSec; 

    double ph1Eta = newMomFst.Eta();
    double ph2Eta = newMomSec.Eta();

    double ph1Phi = newMomFst.Phi();
    double ph2Phi = newMomSec.Phi();
    
    if(print && iVtx == 0 && false ) {
      std::cout<<"     new momenta Et1 = "<<newMomFst.Et()<<std::endl;
      std::cout<<"                 Eta = "<<newMomFst.Eta()<<std::endl; 
      std::cout<<"                 Phi = "<<newMomFst.Phi()<<std::endl; 
      std::cout<<"                  Px = "<<newMomFst.Px()<<std::endl; 
      std::cout<<"                  Py = "<<newMomFst.Py()<<std::endl; 
      std::cout<<"     new momenta Et2 = "<<newMomSec.Et()<<std::endl;
      std::cout<<"                 Eta = "<<newMomSec.Eta()<<std::endl; 
      std::cout<<"                 Phi = "<<newMomSec.Phi()<<std::endl;
      std::cout<<"                  Px = "<<newMomSec.Px()<<std::endl; 
      std::cout<<"                  Py = "<<newMomSec.Py()<<std::endl;  
    }

    FourVectorM totTrkMom;
    for(unsigned int iTrk = 0; iTrk < tVtx->NTracks(); ++iTrk) {
      const Track* tTrk = tVtx->Trk(iTrk);
      //if(tTrk->Pt()<1.) continue;
      double tEta = tTrk->Eta();
      double tPhi = tTrk->Phi();
      double dEta1 = TMath::Abs(tEta-ph1Eta);
      double dEta2 = TMath::Abs(tEta-ph2Eta);
      double dPhi1 = TMath::Abs(tPhi-ph1Phi);
      double dPhi2 = TMath::Abs(tPhi-ph2Phi);
      if(dPhi1 > M_PI) dPhi1 = 2*M_PI - dPhi1;
      if(dPhi2 > M_PI) dPhi2 = 2*M_PI - dPhi2;

      double dR1 = TMath::Sqrt(dEta1*dEta1+dPhi1*dPhi1);
      double dR2 = TMath::Sqrt(dEta2*dEta2+dPhi2*dPhi2);
      
      if( ( iVtx == 0 || iVtx == 1 ) && print) {
        std::cout<<"  Track #"<<iTrk<<std::endl;
        std::cout<<"      pt = "<<tTrk->Pt()<<std::endl;
        std::cout<<"     eta = "<<tTrk->Eta()<<std::endl;
        std::cout<<"     phi = "<<tTrk->Phi()<<std::endl;
        std::cout<<"      px = "<<tTrk->Px()<<std::endl;
        std::cout<<"      py = "<<tTrk->Py()<<std::endl;
        std::cout<<"     dR1 = "<<dR1<<std::endl;
        std::cout<<"     dR2 = "<<dR2<<std::endl;
      } 

      if(dR1 < 0.05 || dR2 < 0.05) continue;

      if(iTrk == 0) totTrkMom = tTrk->Mom4(0);
      else totTrkMom = totTrkMom + tTrk->Mom4(0);

    }

    if(iVtx ==0 && print && false) {
      std::cout<<"   Trk passes cuts: "<<std::endl;
      std::cout<<"            px tot = "<<totTrkMom.Px()<<std::endl;
      std::cout<<"            py tot = "<<totTrkMom.Py()<<std::endl;
    }

    double ptvtx   = totTrkMom.Pt();
    if(iVtx ==0 && print && false)
      std::cout<<"    Total TkPt = "<<ptvtx<<std::endl;
    double pthiggs = higgsMom.Pt();
    if(iVtx ==0 && print && false)
      std::cout<<"    Total H Pt = "<<pthiggs<<std::endl;
    if(iVtx == 0) ptgg = pthiggs;
    ptbal [iVtx]  = (totTrkMom.Px()*(newMomFst.Px()+newMomSec.Px()));
    ptbal [iVtx] += (totTrkMom.Py()*(newMomFst.Py()+newMomSec.Py()));
    ptbal [iVtx]  = -ptbal[iVtx]/pthiggs;
    ptasym[iVtx] = (ptvtx - pthiggs)/(ptvtx + pthiggs);

    if(iVtx ==0 && print && false) {
      std::cout<<"    Results: ptbal = "<<ptbal [iVtx]<<std::endl;
      std::cout<<"            ptasym = "<<ptasym[iVtx]<<std::endl;
    }

    for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
      if(ptbal [iVtx] > ptbal [cVtx])
        ptbal_rank[cVtx]++;
      else
        ptbal_rank[iVtx]++;
      if(ptasym [iVtx] > ptasym [cVtx])
        ptasym_rank[cVtx]++;
      else
        ptasym_rank[iVtx]++;
    }
  }
  

  // compute the total rank
  for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
    if(print) {
      std::cout<<"     Vertex #"<<iVtx<<"  has rank PTBAL "<<ptbal_rank[iVtx]<<"   ("<<ptbal[iVtx]<<")"<<std::endl;
      std::cout<<"     Vertex #"<<iVtx<<"  has rank PTSYM "<<ptasym_rank[iVtx]<<"   ("<<ptasym[iVtx]<<")"<<std::endl;
    }
    ptasym_rank [iVtx] = ptbal_rank [iVtx]*ptasym_rank [iVtx]*(iVtx+1);
    total_rank [iVtx] = 0;
    for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
      if(ptasym_rank [iVtx] > ptasym_rank [cVtx])
        total_rank[iVtx]++;
      else if(ptasym_rank [iVtx] == ptasym_rank [cVtx]) {
        if(ptbal_rank [iVtx] > ptbal_rank [cVtx])
          total_rank[iVtx]++;
        else
          total_rank[cVtx]++;
      }
      else
        total_rank[cVtx]++;
    }
  }

  if(print) std::cout<<std::endl;

  unsigned int bestIdx = 0;
  for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
    if(total_rank[iVtx] == 0) bestIdx = iVtx;
    if(print) {
      std::cout<<"     Vertex #"<<iVtx<<"  has rank "<<total_rank[iVtx]<<std::endl;
    }
  }
  
  delete ptbal_rank  ;
  delete ptasym_rank ;
  delete total_rank  ;
  delete ptbal       ;
  delete ptasym      ;

  //return bestIdx;

  // check if there's a conversion among the pre-selected photons
  const DecayParticle* conv1 = PhotonTools::MatchedCiCConversion(ph1, fConversions);
  const DecayParticle* conv2 = PhotonTools::MatchedCiCConversion(ph2, fConversions);
    
  double zconv  = 0.;
  double dzconv = 0.;
  if(conv1 || conv2) {
    if( conv1 ){
      const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
      zconv  = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
      if( ph1->SCluster()->AbsEta() < 1.5 ) {
        double rho = conv1->Position().Rho();
        if     ( rho < 15. ) dzconv = 0.06;
        else if( rho < 60. ) dzconv = 0.67;
        else                 dzconv = 2.04;
      } else {
        double z = conv1->Position().Z();
        if     ( z < 50. )   dzconv = 0.18;
        else if( z < 100.)   dzconv = 0.61;
        else                 dzconv = 0.99;
      }
    } else if( !conv1 ) {
      const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
      zconv  = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
      if( ph2->SCluster()->AbsEta() < 1.5 ) {
        double rho = conv2->Position().Rho();
        if     ( rho < 15. ) dzconv = 0.06;
        else if( rho < 60. ) dzconv = 0.67;
        else                 dzconv = 2.04;
      } else {
        double z = conv2->Position().Z();
        if     ( z < 50. )   dzconv = 0.18;
        else if( z < 100.)   dzconv = 0.61;
        else                 dzconv = 0.99;
      }
    } else {
      const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
      double z1  = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
      double dz1 = 0.;
      if( ph1->SCluster()->AbsEta() < 1.5 ) {
        double rho = conv1->Position().Rho();
        if     ( rho < 15. ) dz1 = 0.06;
        else if( rho < 60. ) dz1 = 0.67;
        else                 dz1 = 2.04;
      } else {
        double z = conv1->Position().Z();
        if     ( z < 50. )   dz1 = 0.18;
        else if( z < 100.)   dz1 = 0.61;
        else                 dz1 = 0.99;
      }
      const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
      double z2  = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
      double dz2 = 0.;
      if( ph2->SCluster()->AbsEta() < 1.5 ) {
        double rho = conv2->Position().Rho();
        if     ( rho < 15. ) dz2 = 0.06;
        else if( rho < 60. ) dz2 = 0.67;
        else                 dz2 = 2.04;
      } else {
        double z = conv2->Position().Z();
        if     ( z < 50. )   dz2 = 0.18;
        else if( z < 100.)   dz2 = 0.61;
        else                 dz2 = 0.99;
      }
      zconv  = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2 )*(z1/dz1/dz1 + z2/dz2/dz2) ) ;  // weighted average
      dzconv = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2)) ;
    }
  }

      
  // loop over all ranked Vertices and choose the closest to the Conversion one
  int maxVertices = ( ptgg > 30 ? 3 : 5);
  double minDz = -1;    
  for(unsigned int iVtx =0; iVtx < numVertices && (int) iVtx < maxVertices; ++iVtx) {
    if(total_rank[iVtx] < maxVertices) {
      const Vertex* tVtx = fPV->At(iVtx);
      double tDz = TMath::Abs(zconv - tVtx->Z());
      if( (minDz < 0. || tDz < minDz) && tDz < dzconv ) {
        minDz = tDz;
        bestIdx = iVtx;
      }
    }
  }


  return bestIdx;
}

Revision:	1.4
Committed:	Mon Jun 27 12:32:53 2011 UTC (13 years, 10 months ago) by fabstoec
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+477 -200 lines
Log Message:	update Photon CiC selection
#	User	Rev	Content
1	fabstoec	1.1	#include "MitPhysics/Mods/interface/PhotonCiCMod.h"
2			#include "MitAna/DataTree/interface/PhotonCol.h"
3			#include "MitPhysics/Init/interface/ModNames.h"
4			#include "MitPhysics/Utils/interface/IsolationTools.h"
5			#include "MitPhysics/Utils/interface/PhotonTools.h"
6	fabstoec	1.4	#include <TNtuple.h>
7	fabstoec	1.1
8			using namespace mithep;
9
10			ClassImp(mithep::PhotonCiCMod)
11
12			//--------------------------------------------------------------------------------------------------
13			PhotonCiCMod::PhotonCiCMod(const char name, const char title) :
14			BaseMod(name,title),
15			fPhotonBranchName (Names::gkPhotonBrn),
16			fGoodPhotonsName (ModNames::gkGoodPhotonsName),
17			fTrackBranchName (Names::gkTrackBrn),
18			fPileUpDenName (Names::gkPileupEnergyDensityBrn),
19			fElectronName ("Electrons"),
20	bendavid	1.2	fPhotonPtMin(20.0),
21	fabstoec	1.1	fApplySpikeRemoval(kFALSE),
22			fAbsEtaMax(999.99),
23			fPhotons(0),
24			fTracks(0),
25			fPileUpDen(0),
26			fElectrons(0),
27
28	fabstoec	1.4	fPVName(Names::gkPVBeamSpotBrn),
29	fabstoec	1.1	fPV(0),
30	fabstoec	1.4	fPVFromBranch(kTRUE),
31
32			fConversions(0),
33			fConversionName(Names::gkMvfConversionBrn),
34
35			fBeamspot(0)
36	fabstoec	1.1
37			{
38			// Constructor.
39			}
40
41			//--------------------------------------------------------------------------------------------------
42			void PhotonCiCMod::Process()
43			{
44			// Process entries of the tree.
45
46			LoadEventObject(fPhotonBranchName, fPhotons);
47	bendavid	1.2
48			PhotonOArr *GoodPhotons = new PhotonOArr;
49			GoodPhotons->SetName(fGoodPhotonsName);
50	fabstoec	1.4	GoodPhotons->SetOwner(kTRUE);
51
52			Double_t _tRho = -1.;
53			LoadEventObject(fTrackBranchName, fTracks);
54			LoadEventObject(fPileUpDenName, fPileUpDen);
55			LoadEventObject(fElectronName, fElectrons);
56			LoadEventObject(fPVName, fPV);
57			LoadEventObject(fConversionName, fConversions);
58			LoadBranch("BeamSpot");
59	bendavid	1.2
60	fabstoec	1.4	if(fPileUpDen->GetEntries() > 0)
61			_tRho = (Double_t) fPileUpDen->At(0)->Rho();
62
63			_tRho = 0.;
64
65			bool doVtxSelection = true;
66
67			const EventHeader* evtHead = this->GetEventHeader();
68
69			unsigned int evtNum = evtHead->EvtNum();
70			Float_t _evtNum1 = (Float_t) ( (int) (evtNum/10000.) );
71			Float_t _evtNum2 = (Float_t) ( (int) (evtNum % 10000) );
72
73			double evtNumTest = (int) ( ( (double) _evtNum1 )*10000. + (double) _evtNum2 );
74
75			Float_t _runNum = (Float_t) evtHead->RunNum();
76			Float_t _lumiSec = (Float_t) evtHead->LumiSec();
77
78
79			unsigned int numVertices = fPV->GetEntries();
80
81			const BaseVertex bsp = dynamic_cast<const BaseVertex>(fBeamspot->At(0));
82	fabstoec	1.3
83			PhotonOArr* preselPh = new PhotonOArr;
84	fabstoec	1.1
85	fabstoec	1.3	// 1. we do the pre-selection; but keep the non-passing photons in a secont container...
86			for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
87			const Photon *ph = fPhotons->At(i);
88	fabstoec	1.4
89			if(ph->SCluster()->AbsEta()>=2.5 \|\| (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) continue;
90			if(ph->Et() < 20.) continue;
91			if(ph->HadOverEm() > 0.15) continue;
92	fabstoec	1.3	if(ph->AbsEta() < 1.5) {
93	fabstoec	1.4	if(ph->CoviEtaiEta() > 0.013) continue;
94			//if(ph->EcalRecHitIsoDr03() > 10.) continue;
95	fabstoec	1.3	} else {
96	fabstoec	1.4	if(ph->CoviEtaiEta() > 0.03) continue;
97			//if(ph->EcalRecHitIsoDr03() > 10.) continue;
98	fabstoec	1.3	}
99	fabstoec	1.4
100			Bool_t passSpikeRemovalFilter = kTRUE;
101
102			if (ph->SCluster() && ph->SCluster()->Seed()) {
103			if(ph->SCluster()->Seed()->Energy() > 5.0 &&
104			ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
105			passSpikeRemovalFilter = kFALSE;
106			}
107			if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
108
109	fabstoec	1.3	preselPh->Add(ph);
110			}
111	fabstoec	1.1
112	fabstoec	1.4	// sort both by pt... again ;)
113	fabstoec	1.3	preselPh->Sort();
114	fabstoec	1.4
115			unsigned int numPairs = 0;
116			if( preselPh->GetEntries() > 0)
117			numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
118
119			// create all possible pairs of pre-selected photons
120
121			std::vector<unsigned int> idxFst;
122			std::vector<unsigned int> idxSec;
123			std::vector<bool> pairPasses;
124
125			if(numPairs > 0) {
126			for(unsigned int iFst = 0; iFst <preselPh->GetEntries() - 1; ++iFst) {
127			for(unsigned int iSec = iFst + 1; iSec <preselPh->GetEntries(); ++iSec) {
128			idxFst.push_back(iFst);
129			idxSec.push_back(iSec);
130			pairPasses.push_back(true);
131			}
132			}
133			}
134
135			// array to store the index of 'chosen Vtx' for each pair
136			const Vertex** theVtx = new const Vertex*[numPairs];
137
138			// arays to store the Vtx 'fixed' photons
139			Photon** fixPhFst = new Photon*[numPairs];
140			Photon** fixPhSec = new Photon*[numPairs];
141
142			// sotr pair-indices for pairs passing the selection
143			std::vector<unsigned int> passPairs;
144			passPairs.resize(0);
145
146			unsigned int theChosenVtx = 0;
147
148			float kinPh1[20];
149			float kinPh2[20];
150
151			for(int i=0; i<10; ++i) {
152			kinPh1[i] =-99.;
153			kinPh2[i] =-99.;
154			}
155
156			float ptBefore1 = -99.;
157			float ptBefore2 = -99.;
158
159			bool print = false;
160			//if(evtNum == 1677 \|\| evtNum == 1943 \|\| evtNum == 1694) {
161			if(evtNum == 1943 && false) {
162			std::cout<<" ------------------------------------------- "<<std::endl;
163			std::cout<<" printing info for event #"<<evtNum<<std::endl;
164			print = true;
165			}
166
167			for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
168
169			// copy the photons for manipulation
170			fixPhFst[iPair] = new Photon(*preselPh->At(idxFst[iPair]));
171			fixPhSec[iPair] = new Photon(*preselPh->At(idxSec[iPair]));
172
173			// store the vertex for this pair (TODO: conversion Vtx)
174			if(doVtxSelection) {
175			unsigned int iVtx = findBestVertex(fixPhFst[iPair],fixPhSec[iPair],bsp, print);
176			theVtx[iPair] = fPV->At(iVtx);
177			if(iPair == 0) theChosenVtx = iVtx;
178			} else
179			theVtx[iPair] = fPV->At(0);
180
181			if(iPair==0) {
182			ptBefore1 = fixPhFst[iPair]->Pt();
183			ptBefore2 = fixPhSec[iPair]->Pt();
184			}
185
186			// fix the kinematics for both events
187			FourVectorM newMomFst = fixPhFst[iPair]->MomVtx(theVtx[iPair]->Position());
188			FourVectorM newMomSec = fixPhSec[iPair]->MomVtx(theVtx[iPair]->Position());
189			fixPhFst[iPair]->SetMom(newMomFst.X(), newMomFst.Y(), newMomFst.Z(), newMomFst.E());
190			fixPhSec[iPair]->SetMom(newMomSec.X(), newMomSec.Y(), newMomSec.Z(), newMomSec.E());
191
192			if(iPair != 0) {
193			// check if both photons pass the CiC selection
194			bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., false);
195			bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., false);
196			if( pass1 && pass2 )
197			passPairs.push_back(iPair);
198			} else {
199			bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., false, kinPh1);
200			bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., false, kinPh2);
201			if( pass1 && pass2 )
202			passPairs.push_back(iPair);
203			}
204			}
205
206
207			// loop over all passing pairs and find the one with the highest sum Et
208			Photon* phHard = NULL;
209			Photon* phSoft = NULL;
210			double maxSumEt = 0.;
211			for(unsigned int iPair=0; iPair<passPairs.size(); ++iPair){
212			double sumEt = fixPhFst[passPairs[iPair]]->Et();
213			sumEt += fixPhSec[passPairs[iPair]]->Et();
214			if( sumEt > maxSumEt ) {
215			maxSumEt = sumEt;
216			phHard = fixPhFst[passPairs[iPair]];
217			phSoft = fixPhSec[passPairs[iPair]];
218			}
219			}
220
221			if(phHard && phSoft) {
222			GoodPhotons->AddOwned(phHard);
223			GoodPhotons->AddOwned(phSoft);
224			}
225
226			// sort according to pt
227			GoodPhotons->Sort();
228
229			// add to event for other modules to use
230			AddObjThisEvt(GoodPhotons);
231
232			delete preselPh;
233
234
235			bool doFill = (phHard && phSoft);
236			Float_t _mass = ( doFill ? (phHard->Mom()+phSoft->Mom()).M() : -100.);
237			Float_t _ptgg = ( doFill ? (phHard->Mom()+phSoft->Mom()).Pt() : -100.);
238
239
240			Float_t fillEvent[] = { _mass,
241			_ptgg,
242			_evtNum1,
243			_evtNum2,
244			_runNum,
245			_lumiSec,
246			(float) theChosenVtx,
247			(float) numPairs,
248			kinPh1[0],
249			kinPh1[1],
250			kinPh1[2],
251			kinPh1[3],
252			kinPh1[4],
253			kinPh1[5],
254			kinPh1[6],
255			kinPh1[7],
256			kinPh1[8],
257			kinPh1[9],
258			kinPh1[10],
259			kinPh1[11],
260			kinPh1[12],
261			kinPh1[13],
262			kinPh1[14],
263			kinPh1[15],
264			kinPh1[16],
265			kinPh1[17],
266			kinPh1[18],
267			kinPh1[19],
268			kinPh2[0],
269			kinPh2[1],
270			kinPh2[2],
271			kinPh2[3],
272			kinPh2[4],
273			kinPh2[5],
274			kinPh2[6],
275			kinPh2[7],
276			kinPh2[8],
277			kinPh2[9],
278			kinPh2[10],
279			kinPh2[11],
280			kinPh2[12],
281			kinPh2[13],
282			kinPh2[14],
283			kinPh2[15],
284			kinPh2[16],
285			kinPh2[17],
286			kinPh2[18],
287			kinPh2[19],
288			ptBefore1,
289			ptBefore2
290			};
291
292			hCiCTuple->Fill(fillEvent);
293
294			return;
295
296			}
297
298			//--------------------------------------------------------------------------------------------------
299			void PhotonCiCMod::SlaveBegin()
300			{
301			// Run startup code on the computer (slave) doing the actual analysis. Here,
302			// we just request the photon collection branch.
303
304			ReqEventObject(fPhotonBranchName, fPhotons, kTRUE);
305			ReqEventObject(fTrackBranchName, fTracks, kTRUE);
306			ReqEventObject(fElectronName, fElectrons, kTRUE);
307			ReqEventObject(fPileUpDenName, fPileUpDen, kTRUE);
308			ReqEventObject(fPVName, fPV, fPVFromBranch);
309			ReqEventObject(fConversionName, fConversions,kTRUE);
310			ReqBranch("BeamSpot",fBeamspot);
311
312
313			hCiCTuple = new TNtuple("hCiCTuple","hCiCTuple","mass:ptgg:evtnum1:evtnum2:runnum:lumisec:ivtx:npairs:ph1Iso1:ph1Iso2:ph1Iso3:ph1Cov:ph1HoE:ph1R9:ph1DR:ph1Pt:ph1Eta:ph1Phi:ph1Eiso3:ph1Eiso4:ph1Hiso4:ph1TisoA:ph1TisoW:ph1Tiso:ph1Et:ph1E:ph1Pass:ph1Cat:ph2Iso1:ph2Iso2:ph2Iso3:ph2Cov:ph2HoE:ph2R9:ph2DR:ph2Pt:ph2Eta:ph2Phi:ph2Eiso3:ph2Eiso4:ph2Hiso4:ph2TisoA:ph2TisoW:ph2Tiso:ph2Et:ph2E:ph2Pass:ph2Cat:ph1UPt:ph2UPt");
314	fabstoec	1.3
315	fabstoec	1.4	AddOutput(hCiCTuple);
316
317			}
318
319			// return the index of the bext vertex
320			unsigned int PhotonCiCMod::findBestVertex(Photon* ph1, Photon* ph2, const BaseVertex *bsp, bool print) {
321	fabstoec	1.3
322			// loop over all vertices and assigne the ranks
323			int* ptbal_rank = new int[fPV->GetEntries()];
324			int* ptasym_rank = new int[fPV->GetEntries()];
325			int* total_rank = new int[fPV->GetEntries()];
326			double* ptbal = new double[fPV->GetEntries()];
327			double* ptasym = new double[fPV->GetEntries()];
328	fabstoec	1.4
329			unsigned int numVertices = fPV->GetEntries();
330
331			double ptgg = 0.; // stored for later in the conversion
332	fabstoec	1.3
333	fabstoec	1.4	if(print) {
334			std::cout<<" --------------------------------- "<<std::endl;
335			std::cout<<" looking for Vtx for photon pair "<<std::endl;
336			std::cout<<" pt 1 = "<<ph1->Et()<<std::endl;
337			std::cout<<" pt 2 = "<<ph2->Et()<<std::endl;
338			std::cout<<" among #"<<numVertices<<" Vertices."<<std::endl;
339			}
340
341
342			for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
343			if(print)
344			std::cout<<std::endl<<" Vertex #"<<iVtx<<std::endl;
345
346			const Vertex* tVtx = fPV->At(iVtx);
347	fabstoec	1.3	ptbal [iVtx] = 0.0;
348			ptasym[iVtx] = 0.0;
349			ptbal_rank [iVtx] = 1;
350			ptasym_rank[iVtx] = 1;
351	fabstoec	1.4
352			// compute the photon momenta with respect to this Vtx
353			FourVectorM newMomFst = ph1->MomVtx(tVtx->Position());
354			FourVectorM newMomSec = ph2->MomVtx(tVtx->Position());
355
356			FourVectorM higgsMom = newMomFst+newMomSec;
357
358			double ph1Eta = newMomFst.Eta();
359			double ph2Eta = newMomSec.Eta();
360
361			double ph1Phi = newMomFst.Phi();
362			double ph2Phi = newMomSec.Phi();
363
364			if(print && iVtx == 0 && false ) {
365			std::cout<<" new momenta Et1 = "<<newMomFst.Et()<<std::endl;
366			std::cout<<" Eta = "<<newMomFst.Eta()<<std::endl;
367			std::cout<<" Phi = "<<newMomFst.Phi()<<std::endl;
368			std::cout<<" Px = "<<newMomFst.Px()<<std::endl;
369			std::cout<<" Py = "<<newMomFst.Py()<<std::endl;
370			std::cout<<" new momenta Et2 = "<<newMomSec.Et()<<std::endl;
371			std::cout<<" Eta = "<<newMomSec.Eta()<<std::endl;
372			std::cout<<" Phi = "<<newMomSec.Phi()<<std::endl;
373			std::cout<<" Px = "<<newMomSec.Px()<<std::endl;
374			std::cout<<" Py = "<<newMomSec.Py()<<std::endl;
375			}
376
377			FourVectorM totTrkMom;
378			for(unsigned int iTrk = 0; iTrk < tVtx->NTracks(); ++iTrk) {
379			const Track* tTrk = tVtx->Trk(iTrk);
380			//if(tTrk->Pt()<1.) continue;
381			double tEta = tTrk->Eta();
382			double tPhi = tTrk->Phi();
383			double dEta1 = TMath::Abs(tEta-ph1Eta);
384			double dEta2 = TMath::Abs(tEta-ph2Eta);
385			double dPhi1 = TMath::Abs(tPhi-ph1Phi);
386			double dPhi2 = TMath::Abs(tPhi-ph2Phi);
387			if(dPhi1 > M_PI) dPhi1 = 2*M_PI - dPhi1;
388			if(dPhi2 > M_PI) dPhi2 = 2*M_PI - dPhi2;
389
390			double dR1 = TMath::Sqrt(dEta1dEta1+dPhi1dPhi1);
391			double dR2 = TMath::Sqrt(dEta2dEta2+dPhi2dPhi2);
392
393			if( ( iVtx == 0 \|\| iVtx == 1 ) && print) {
394			std::cout<<" Track #"<<iTrk<<std::endl;
395			std::cout<<" pt = "<<tTrk->Pt()<<std::endl;
396			std::cout<<" eta = "<<tTrk->Eta()<<std::endl;
397			std::cout<<" phi = "<<tTrk->Phi()<<std::endl;
398			std::cout<<" px = "<<tTrk->Px()<<std::endl;
399			std::cout<<" py = "<<tTrk->Py()<<std::endl;
400			std::cout<<" dR1 = "<<dR1<<std::endl;
401			std::cout<<" dR2 = "<<dR2<<std::endl;
402			}
403
404			if(dR1 < 0.05 \|\| dR2 < 0.05) continue;
405
406			if(iTrk == 0) totTrkMom = tTrk->Mom4(0);
407			else totTrkMom = totTrkMom + tTrk->Mom4(0);
408
409			}
410
411			if(iVtx ==0 && print && false) {
412			std::cout<<" Trk passes cuts: "<<std::endl;
413			std::cout<<" px tot = "<<totTrkMom.Px()<<std::endl;
414			std::cout<<" py tot = "<<totTrkMom.Py()<<std::endl;
415			}
416
417			double ptvtx = totTrkMom.Pt();
418			if(iVtx ==0 && print && false)
419			std::cout<<" Total TkPt = "<<ptvtx<<std::endl;
420			double pthiggs = higgsMom.Pt();
421			if(iVtx ==0 && print && false)
422			std::cout<<" Total H Pt = "<<pthiggs<<std::endl;
423			if(iVtx == 0) ptgg = pthiggs;
424			ptbal [iVtx] = (totTrkMom.Px()*(newMomFst.Px()+newMomSec.Px()));
425			ptbal [iVtx] += (totTrkMom.Py()*(newMomFst.Py()+newMomSec.Py()));
426			ptbal [iVtx] = -ptbal[iVtx]/pthiggs;
427			ptasym[iVtx] = (ptvtx - pthiggs)/(ptvtx + pthiggs);
428
429			if(iVtx ==0 && print && false) {
430			std::cout<<" Results: ptbal = "<<ptbal [iVtx]<<std::endl;
431			std::cout<<" ptasym = "<<ptasym[iVtx]<<std::endl;
432			}
433
434			for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
435	fabstoec	1.3	if(ptbal [iVtx] > ptbal [cVtx])
436	fabstoec	1.4	ptbal_rank[cVtx]++;
437	fabstoec	1.3	else
438	fabstoec	1.4	ptbal_rank[iVtx]++;
439	fabstoec	1.3	if(ptasym [iVtx] > ptasym [cVtx])
440	fabstoec	1.4	ptasym_rank[cVtx]++;
441	fabstoec	1.3	else
442	fabstoec	1.4	ptasym_rank[iVtx]++;
443	fabstoec	1.3	}
444			}
445	fabstoec	1.4
446	fabstoec	1.1
447	fabstoec	1.3	// compute the total rank
448	fabstoec	1.4	for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
449			if(print) {
450			std::cout<<" Vertex #"<<iVtx<<" has rank PTBAL "<<ptbal_rank[iVtx]<<" ("<<ptbal[iVtx]<<")"<<std::endl;
451			std::cout<<" Vertex #"<<iVtx<<" has rank PTSYM "<<ptasym_rank[iVtx]<<" ("<<ptasym[iVtx]<<")"<<std::endl;
452			}
453			ptasym_rank [iVtx] = ptbal_rank [iVtx]ptasym_rank [iVtx](iVtx+1);
454	fabstoec	1.3	total_rank [iVtx] = 0;
455	fabstoec	1.4	for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
456			if(ptasym_rank [iVtx] > ptasym_rank [cVtx])
457	fabstoec	1.3	total_rank[iVtx]++;
458	fabstoec	1.4	else if(ptasym_rank [iVtx] == ptasym_rank [cVtx]) {
459			if(ptbal_rank [iVtx] > ptbal_rank [cVtx])
460			total_rank[iVtx]++;
461			else
462			total_rank[cVtx]++;
463			}
464	fabstoec	1.3	else
465			total_rank[cVtx]++;
466			}
467	bendavid	1.2	}
468	fabstoec	1.1
469	fabstoec	1.4	if(print) std::cout<<std::endl;
470
471			unsigned int bestIdx = 0;
472			for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
473			if(total_rank[iVtx] == 0) bestIdx = iVtx;
474			if(print) {
475			std::cout<<" Vertex #"<<iVtx<<" has rank "<<total_rank[iVtx]<<std::endl;
476			}
477			}
478	fabstoec	1.3
479			delete ptbal_rank ;
480			delete ptasym_rank ;
481			delete total_rank ;
482			delete ptbal ;
483			delete ptasym ;
484	fabstoec	1.1
485	fabstoec	1.4	//return bestIdx;
486
487			// check if there's a conversion among the pre-selected photons
488			const DecayParticle* conv1 = PhotonTools::MatchedCiCConversion(ph1, fConversions);
489			const DecayParticle* conv2 = PhotonTools::MatchedCiCConversion(ph2, fConversions);
490
491			double zconv = 0.;
492			double dzconv = 0.;
493			if(conv1 \|\| conv2) {
494			if( conv1 ){
495			const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
496			zconv = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
497			if( ph1->SCluster()->AbsEta() < 1.5 ) {
498			double rho = conv1->Position().Rho();
499			if ( rho < 15. ) dzconv = 0.06;
500			else if( rho < 60. ) dzconv = 0.67;
501			else dzconv = 2.04;
502			} else {
503			double z = conv1->Position().Z();
504			if ( z < 50. ) dzconv = 0.18;
505			else if( z < 100.) dzconv = 0.61;
506			else dzconv = 0.99;
507			}
508			} else if( !conv1 ) {
509			const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
510			zconv = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
511			if( ph2->SCluster()->AbsEta() < 1.5 ) {
512			double rho = conv2->Position().Rho();
513			if ( rho < 15. ) dzconv = 0.06;
514			else if( rho < 60. ) dzconv = 0.67;
515			else dzconv = 2.04;
516			} else {
517			double z = conv2->Position().Z();
518			if ( z < 50. ) dzconv = 0.18;
519			else if( z < 100.) dzconv = 0.61;
520			else dzconv = 0.99;
521			}
522			} else {
523			const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
524			double z1 = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
525			double dz1 = 0.;
526			if( ph1->SCluster()->AbsEta() < 1.5 ) {
527			double rho = conv1->Position().Rho();
528			if ( rho < 15. ) dz1 = 0.06;
529			else if( rho < 60. ) dz1 = 0.67;
530			else dz1 = 2.04;
531			} else {
532			double z = conv1->Position().Z();
533			if ( z < 50. ) dz1 = 0.18;
534			else if( z < 100.) dz1 = 0.61;
535			else dz1 = 0.99;
536			}
537			const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
538			double z2 = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
539			double dz2 = 0.;
540			if( ph2->SCluster()->AbsEta() < 1.5 ) {
541			double rho = conv2->Position().Rho();
542			if ( rho < 15. ) dz2 = 0.06;
543			else if( rho < 60. ) dz2 = 0.67;
544			else dz2 = 2.04;
545			} else {
546			double z = conv2->Position().Z();
547			if ( z < 50. ) dz2 = 0.18;
548			else if( z < 100.) dz2 = 0.61;
549			else dz2 = 0.99;
550			}
551			zconv = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2 )*(z1/dz1/dz1 + z2/dz2/dz2) ) ; // weighted average
552			dzconv = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2)) ;
553	fabstoec	1.1	}
554	fabstoec	1.4	}
555	fabstoec	1.1
556	fabstoec	1.4
557			// loop over all ranked Vertices and choose the closest to the Conversion one
558			int maxVertices = ( ptgg > 30 ? 3 : 5);
559			double minDz = -1;
560			for(unsigned int iVtx =0; iVtx < numVertices && (int) iVtx < maxVertices; ++iVtx) {
561			if(total_rank[iVtx] < maxVertices) {
562			const Vertex* tVtx = fPV->At(iVtx);
563			double tDz = TMath::Abs(zconv - tVtx->Z());
564			if( (minDz < 0. \|\| tDz < minDz) && tDz < dzconv ) {
565			minDz = tDz;
566			bestIdx = iVtx;
567			}
568			}
569	fabstoec	1.1	}
570
571
572
573	fabstoec	1.4	return bestIdx;
574	fabstoec	1.1	}
575	fabstoec	1.3