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),

  fDataEnCorr_EB_hR9(0.0049),
  fDataEnCorr_EB_lR9(0.0011),
  fDataEnCorr_EE_hR9(0.0057),
  fDataEnCorr_EE_lR9(0.0039)

{
  // 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)->RhoRandomLowEta();

  //_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]));

    // if this is Data, scale the energy...

    // 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[] = { _tRho,
                          _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","rho: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);
  if( conv1 && conv1->Prob() < 0.0005) conv1 = NULL;
  if( conv2 && conv2->Prob() < 0.0005) conv2 = NULL;

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