Mods/src/ConversionRemoval.cc

// $Id: ConversionRemoval.cc,v 1.3 2009/01/12 10:26:47 bendavid Exp $
#include <TH1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TCanvas.h>

#include "MitConversions/Mods/interface/ConversionRemoval.h"
#include "MitAna/DataTree/interface/Electron.h"
#include "MitAna/DataTree/interface/Track.h"
#include "MitAna/DataCont/interface/ObjArray.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitCommon/MathTools/interface/MathUtils.h"

using namespace mithep;

ClassImp(mithep::ConversionRemoval)

//__________________________________________________________________________________________________
ConversionRemoval::ConversionRemoval(const char *name, const char *title) :
  BaseMod    (name,title),
  fPtMin   (0),
  fEtaMax  (0),
  fProbMin (0),
  fNHitsMin(0),
  fNStereoHitsMin(0),
  fNPixelHitsMin(0),
  fExcludePXB1(false),
  fMassMax(0),
  fMissedHitsMax(0),
  fWrongHitsMax(0),
  fElectronPtMin(0),
  fElectronEtaMax(0),
  fComputeEff(kFALSE),
  fHitBasedMatching(kFALSE),
  fRequireCleanElectron(kFALSE),
  fTrackQualityFirstOnly(kFALSE),
  fMatchWe(kFALSE),
  fMatchSimConv(kFALSE),
  fRequireGoodPv(kFALSE),
  fMaxTrackDzPv(999.99),
  fMinClusterVtxQual(0.0),
  fMaxNSharedHits(999),
  fMinNSharedHits(0),
  fTracks     (0),
  fGsfTracks  (0),
  fInOutTracks (0),
  fOutInTracks (0),
  fElectrons  (0),
  fTrkElectrons(0),
  fPhotons (0),
  fMCParticles (0),
  fMvfConversions (0),
  fMvfConversionsUnconstrained (0),
  fPrimaryVertexes(0),
  fEvtSel(0),
  fCleanElectrons(0),
  fMCPartName (Names::gkMCPartBrn),
  fTrackName  (Names::gkTrackBrn),
  fConvElectronName   ("Electrons"),
  fPhotonName   (Names::gkPhotonBrn),
  fMvfConvName("ConversionRemoval"),
  fMvfConvUnconstrainedName("ConversionRemovalUnconstrained"),
  fElectronName(Names::gkElectronBrn),
  fTrkElectronName("TrackerElectrons"),
  fGoodElectronsName("ConversionElectrons"),
  fGoodTwoClusterElectronsName("TwoClusterConversionElectrons"),
  fBadElectronsName("ConversionRemovedElectrons"),
  fFakingConversionsName("fakingConversions"),
  fRemovedElectronsName("removedElectrons"),
  fAllConversionElectronsName("allConversionElectrons"),
  hConversionRadius   (0),
  hConversionRadiusBarrel(0),
  hConversionZ(0),
  hConversionRPhi (0),
  hConversionRZ(0),
  hConversionFoldedRZ(0),
  hGammaPt (0),
  hGammaEta (0),
  hGammaPhi(0),
  hGammaMass (0),
  hSimMatchedGammaMass(0),
  hConvProb (0),
  hConvChi2 (0),
  hSimMatchedConvChi2(0),
  hConvDCotTheta (0),
  hConvDCotThetaPreCut (0),
  hConvEOverP (0),
  hElectronPt (0),
  hMinElectronPt (0),
  hSimMinElectronPt (0),
  hSimMatchedMinElectronPt (0),
  hSimMatchedMinElectronSimPt (0),
  hElectronEta (0),
  hEPairMass (0),
  hEPairDeltaPhi (0),
  hGenNumDaughters (0),
  hGenDaughterPt (0),
  hGenDecayRadius (0),
  hSimPt (0),
  hSimEta(0),
  hSimConvRadius (0),
  hSimConvPosition(0),
  hTrackSimMatchType (0),
  hParentSimMatchType (0),
  hSimMatchedConvRadius (0),
  hSimMatchedConvProb (0),
  hSimMatchedEPairMass (0),
  hSimMatchedEPairDeltaPhi (0),
  hSimMatchedSimConvRadius (0),
  hSimMatchedConvResolution (0),
  hSimMatchedConvPhiRes (0),
  hSimMatchedConvEOverP (0),
  hSimMatchedElectronPt (0),
  hSimMatchedGammaPt (0),
  hSimMatchedGammaEta (0),
  hSimMatchedSimGammaPt (0),
  hSimMatchedSimGammaEta (0),
  hTrackChi2          (0),
  hTrackProb          (0),
  hTrackNHits         (0),
  hTrackNWrongHits(0),
  hTrackNMissingHits(0),
  hNSharedHits(0),
  hTrackNHitsProb     (0),
  hTrackD0 (0),
  hTrackDzPv(0),
  hUnMatchedTrackChi2          (0),
  hUnMatchedTrackProb          (0),
  hUnMatchedTrackNHits         (0),
  hUnMatchedTrackNHitsProb     (0),
  hUnMatchedTrackD0 (0),
  hSimMatchedTrackChi2          (0),
  hSimMatchedTrackProb          (0),
  hSimMatchedTrackNHits         (0),
  hSimMatchedTrackNHitsProb     (0),
  hSimMatchedTrackD0 (0),
  hLxy(0),
  hLxyOverLxyErr(0),
  hLz(0),
  hLzOverLzErr(0),
  hDxy(0),
  hSimMatchedLxy(0),
  hSimMatchedLz(0),
  hSimMatchedDxy(0),
  hNTracks(0),
  hSimMatchedNTracks(0),
  hIsolation(0),
  hSimMatchedIsolation(0),
  hNConversions(0),
  hDoubleConvMass(0),
  hConvGammaMass(0),
  hPiPiMass(0),
  hAllConversionElectronChargeAssign(0),
  hRemovedElectronChargeAssign(0),
  hTrackAlgo(0),
  hCharge(0),
  hPt1Pt2(0),
  hClusVtxDiff(0),
  hClusVtxQual(0)
{
  // Constructor.
}

//__________________________________________________________________________________________________
void ConversionRemoval::Begin()
{
  // Run startup code on the client machine. For this module, we dont do anything here.
}

//__________________________________________________________________________________________________
void ConversionRemoval::SlaveBegin()
{
  // Run startup code on the computer (slave) doing the actual analysis. Here, we typically
  // initialize histograms and other analysis objects and request branches. For this module, we
  // request a branch of the MitTree.

  // Request the branches
  ReqBranch(fTrackName,  fTracks);
  ReqBranch("ConversionInOutTracks", fInOutTracks);
  ReqBranch("ConversionOutInTracks", fOutInTracks);
  //ReqBranch(fElectronName,   fElectrons);
  //ReqBranch(fTrkElectronName,   fTrkElectrons);
  ReqBranch(fPhotonName,   fPhotons);
  ReqBranch(fMCPartName,   fMCParticles);
  ReqBranch(fMvfConvName, fMvfConversions);
  ReqBranch(fMvfConvUnconstrainedName, fMvfConversionsUnconstrained);
  ReqBranch("PrimaryVertexes",fPrimaryVertexes);
  ReqBranch("EvtSelData",fEvtSel);

  // Book our histograms
  
  hConversionRPhi = new TH2F("Photon Conversion Position R-Phi","Conversion Position R-Phi (cm)", 100, -100.0,100.0,100,-100.0,100.0);
  AddOutput(hConversionRPhi);
  
  hConversionRZ = new TH2F("Photon Conversion Position R-Z","Conversion Position R-Z (cm)", 1000, -100.0,100.0,1000,-100.0,100.0);
  AddOutput(hConversionRZ);
  
  hConversionFoldedRZ = new TH2F("hConversionFoldedRZ","Conversion Position R-Z (cm)", 1000, -100.0,100.0,500,-0.0,100.0);
  AddOutput(hConversionFoldedRZ);
  
  hSimConvPosition = new TH2F("Sim Photon Conversion Position","Sim Photon Conversion Position", 200, -200.0,200.0,200,-200.0,200.0);
  AddOutput(hSimConvPosition);
  
  
  //TH1::SetDefaultSumw2(kTRUE);
  
  hConversionRadius = new TH1F("Photon Conversion Radius","Radius (cm)",800,0,200.0);
  AddOutput(hConversionRadius);
  
  hConversionRadiusBarrel = new TH1F("hConversionRadiusBarrel","Radius (cm)",800,0,200.0);
  AddOutput(hConversionRadiusBarrel);
  
  hConversionZ = new TH1F("Photon Conversion Z position", "Conversion Z position (cm)", 200, -200.0, 200.0);
  AddOutput(hConversionZ);
  
  hGammaPt = new TH1F("Photon Pt","Photon Pt (GeV)",400,0,200.0);
  AddOutput(hGammaPt);
  
  hGammaEta = new TH1F("Photon Eta","Photon Eta",200,-5.0,5.0);
  AddOutput(hGammaEta);
  
  hGammaPhi = new TH1F("hGammaPhi","Photon Phi",200,-1.0*TMath::Pi(),1.0*TMath::Pi());
  AddOutput(hGammaPhi);
  
  
  hGammaMass = new TH1F("Photon Invariant Mass","Photon Invariant Mass",200,0,2.0);
  AddOutput(hGammaMass);

  hSimMatchedGammaMass = new TH1F("SimMatched Photon Invariant Mass","SimMatched Photon Invariant Mass",200,0,20.0);
  AddOutput(hSimMatchedGammaMass);  

  hConvProb = new TH1F("Conversion Vertex fit probability", "Vertex Fit Probability",200,0.0,1.0);
  AddOutput(hConvProb);
  
  hConvChi2 = new TH1F("Conversion Vertex fit Chi squared", "Vertex Fit Chi Squared",2000,0.0,100.0);
  AddOutput(hConvChi2);

  hSimMatchedConvChi2 = new TH1F("Sim Matched Conversion Vertex fit Chi squared", "Sim Matched Vertex Fit Chi Squared",2000,0.0,100.0);
  AddOutput(hSimMatchedConvChi2);
  
  hConvDCotTheta = new TH1F("Conversion DCotTheta", "DCotTheta",10001,-10.0,10.0);
  AddOutput(hConvDCotTheta);
  
  hConvDCotThetaPreCut = new TH1F("Conversion DCotTheta before fit prob cut", "DCotThetaPreCut",1000,-10.0,10.0);
  AddOutput(hConvDCotThetaPreCut);
  
  hConvEOverP = new TH1F("Conversion E over P", "Conversion E over P",200,0.0,50.0);
  AddOutput(hConvEOverP);
  
  hElectronPt = new TH1F("Electron Pt","Electron Pt (GeV)",400,0,200.0);
  AddOutput(hElectronPt);
  
  hMinElectronPt = new TH1F("Min Electron Pt","Min Electron Pt (GeV)",2000,0,200.0);
  AddOutput(hMinElectronPt);
  
  hSimMatchedMinElectronPt = new TH1F("Sim Matched Min Electron Pt","Sim Matched Min Electron Pt (GeV)",2000,0,200.0);
  AddOutput(hSimMatchedMinElectronPt);
  
  hSimMatchedMinElectronSimPt = new TH1F("Sim Matched Min Electron Sim Pt","Sim Matched Min Electron Sim Pt (GeV)",2000,0,200.0);
  AddOutput(hSimMatchedMinElectronSimPt);
  
  hSimMinElectronPt = new TH1F("Sim Min Electron Pt","Sim Min Electron Pt (GeV)",2000,0,200.0);
  AddOutput(hSimMinElectronPt);
  
  hElectronEta = new TH1F("Electron Eta","Electron Eta",200,-5.0,5.0);
  AddOutput(hElectronEta);
  
  hEPairMass = new TH1F("Electron Pair Invariant Mass","Electron Pair Invariant Mass",200,0.0,10);
  AddOutput(hEPairMass);
  
  hEPairDeltaPhi = new TH1F("Electron Pair Delta Phi","Electron Pair Delta Phi",200,-4.0,4.0);
  AddOutput(hEPairDeltaPhi);
  
  hGenNumDaughters = new TH1F("Number of Decay daughters - gen level","Number of Decay daughters - gen level",10,-0.5,9.5);
  AddOutput(hGenNumDaughters);
  
  hGenDaughterPt = new TH1F("Pt of decay daugthers - gen level","Pt of decay daugthers - gen level",2000,0,200.0);
  AddOutput(hGenDaughterPt);
  
  hGenDecayRadius = new TH1F("Decay vertex radius - gen level","Decay vertex radius - gen level",800,0,200.0);
  AddOutput(hGenDecayRadius);
  
  hSimPt = new TH1F("Pt - sim level","Pt - sim level",2000,0,200.0);
  AddOutput(hSimPt);

  hSimEta = new TH1F("Sim Photon Eta","Photon Eta",200,-5.0,5.0);
  AddOutput(hSimEta);
  
  hSimConvRadius = new TH1F("Sim Photon Conversion Radius","Sim Photon Conversion Radius",800,0,200.0);
  AddOutput(hSimConvRadius);
  
//   hTrackD0 = new TH1F("general track D0","general track D0",2000,0,200.0);
//   AddOutput(hTrackD0);
  
  hTrackSimMatchType = new TH1F("Conversion Electron matched sim pdg","Conversion Electron matched sim pdg",10001,-5000.5,5000.5);
  AddOutput(hTrackSimMatchType);
  
  hParentSimMatchType = new TH1F("Conversion Parent matched sim pdg","Conversion Parent matched sim pdg",10001,-5000.5,5000.5);
  AddOutput(hParentSimMatchType);
  
  hSimMatchedConvRadius = new TH1F("Sim Matched Photon Conversion Radius","Sim Matched Photon Conversion Radius",800,0,200.0);
  AddOutput(hSimMatchedConvRadius);
  
  hSimMatchedConvProb = new TH1F("Sim Matched Conversion Vertex fit probability", "Sim Matched Vertex Fit Probability",200,0.0,1.0);
  AddOutput(hSimMatchedConvProb);
  
  hSimMatchedEPairMass = new TH1F("Sim Matched Electron Pair Invariant Mass","Sim Matched Electron Pair Invariant Mass",200,0.0,10);
  AddOutput(hSimMatchedEPairMass);
  
  hSimMatchedEPairDeltaPhi = new TH1F("Sim Matched Electron Pair Delta Phi","Sim Matched Electron Pair Delta Phi",200,-4.0,4.0);
  AddOutput(hSimMatchedEPairDeltaPhi);
  
  hSimMatchedSimConvRadius = new TH1F("Sim Matched Photon Sim Conversion Radius","Sim Matched Photon Conversion Sim Radius",800,0,200.0);
  AddOutput(hSimMatchedSimConvRadius);
  
  hSimMatchedConvResolution = new TH1F("Sim Matched Photon Conversion Radius Resolution","Sim Matched Photon Conversion Radius Resolution",300,-30.0,30.0);
  AddOutput(hSimMatchedConvResolution);
  
  hSimMatchedConvPhiRes = new TH1F("Sim Matched Photon Conversion Phi Resolution","Sim Matched Photon Conversion Phi Resolution",1000,-3.5,3.5);
  AddOutput(hSimMatchedConvPhiRes);
  
  hSimMatchedConvEOverP = new TH1F("Sim Matched Conversion E over P", "Sim Matched Conversion E over P",200,0.0,50.0);
  AddOutput(hSimMatchedConvEOverP);
  
  hSimMatchedElectronPt = new TH1F("Sim Matched Electron Pt","Sim Matched Electron Pt (GeV)",2000,0,200.0);
  AddOutput(hSimMatchedElectronPt);
  
  hSimMatchedGammaPt = new TH1F("Sim Matched Photon Pt","Sim Matched Photon Pt (GeV)",2000,0,200.0);
  AddOutput(hSimMatchedGammaPt);
  
  hSimMatchedGammaEta = new TH1F("Sim Matched Photon Eta","Sim Matched Photon Eta",200,-5.0,5.0);
  AddOutput(hSimMatchedGammaEta);
  
  hSimMatchedSimGammaPt = new TH1F("Sim Matched Sim Photon Pt","Sim Matched Sim Photon Pt (GeV)",2000,0,200.0);
  AddOutput(hSimMatchedSimGammaPt);
  
  hSimMatchedSimGammaEta = new TH1F("Sim Matched Sim Photon Eta","Sim Matched Sim Photon Eta",200,-5.0,5.0);
  AddOutput(hSimMatchedSimGammaEta);
  
  hTrackChi2 = new TH1F("Conversion Track RChi squared", "Track RChi Squared",200,0.0,10.0);
  AddOutput(hTrackChi2);
  
  hTrackProb = new TH1F("Conversion Track fit probability", "Track Fit Probability",200,0.0,1.0);
  AddOutput(hTrackProb);
  
  hTrackNHits = new TH1F("Conversion Track NHits", "Track NHits",37,-0.5,36.5);
  AddOutput(hTrackNHits);
  
  hTrackNMissingHits = new TH1F("hTrackNMissingHits", "Track NMissingHits",37,-0.5,36.5);
  AddOutput(hTrackNMissingHits);
  
  hTrackNWrongHits = new TH1F("hTrackNWrongHits", "Track NWrongHits",37,-0.5,36.5);
  AddOutput(hTrackNWrongHits);
    
  hNSharedHits = new TH1F("hNSharedHits", "hNSharedHits",37,-0.5,36.5);
  AddOutput(hNSharedHits);
  
  hTrackNHitsProb = new TH2F("Track NHits-RChi2","Track NHits-RChi2", 37,-0.5,36.5,50,0.0,10.0);
  AddOutput(hTrackNHitsProb);
  
  hTrackD0 = new TH1F("Conversion track D0","Conversion track D0",100,-5.0,5.0);
  AddOutput(hTrackD0);

  hTrackDzPv = new TH1F("hTrackDzPv","Conversion track DZ wrt PV",100,-20.0,20.0);
  AddOutput(hTrackDzPv);
  
  hUnMatchedTrackChi2 = new TH1F("UnMatched Conversion Track RChi squared", "UnMatched Track RChi Squared",200,0.0,10.0);
  AddOutput(hUnMatchedTrackChi2);
  
  hUnMatchedTrackProb = new TH1F("UnMatched Conversion Track fit probability", "UnMatched Track Fit Probability",200,0.0,1.0);
  AddOutput(hUnMatchedTrackProb);
  
  hUnMatchedTrackNHits = new TH1F("UnMatched Conversion Track NHits", "UnMatched Track NHits",37,-0.5,36.5);
  AddOutput(hUnMatchedTrackNHits);
  
  hUnMatchedTrackNHitsProb = new TH2F("UnMatched Track NHits-RChi2","UnMatched Track NHits-RChi2", 37,-0.5,36.5,50,0.0,10.0);
  AddOutput(hUnMatchedTrackNHitsProb);
  
  hUnMatchedTrackD0 = new TH1F("UnMatched Conversion track D0","UnMatched Conversion track D0",100,-5.0,5.0);
  AddOutput(hUnMatchedTrackD0);
  
  hSimMatchedTrackChi2 = new TH1F("SimMatched Conversion Track RChi squared", "SimMatched Track RChi Squared",200,0.0,10.0);
  AddOutput(hSimMatchedTrackChi2);
  
  hSimMatchedTrackProb = new TH1F("SimMatched Conversion Track fit probability", "SimMatched Track Fit Probability",200,0.0,1.0);
  AddOutput(hSimMatchedTrackProb);
  
  hSimMatchedTrackNHits = new TH1F("SimMatched Conversion Track NHits", "SimMatched Track NHits",37,-0.5,36.5);
  AddOutput(hSimMatchedTrackNHits);
  
  hSimMatchedTrackNHitsProb = new TH2F("SimMatched Track NHits-RChi2","SimMatched Track NHits-RChi2", 37,-0.5,36.5,50,0.0,10.0);
  AddOutput(hSimMatchedTrackNHitsProb);
  
  hSimMatchedTrackD0 = new TH1F("SimMatched Conversion track D0","SimMatched Conversion track D0",100,-5.0,5.0);
  AddOutput(hSimMatchedTrackD0);
  
  hLxy = new TH1F("Conversion Lxy","Conversion Lxy",240,-30.0,30.0);
  AddOutput(hLxy);
  
  hLxyOverLxyErr = new TH1F("Conversion Lxy/LxyErr","Conversion Lxy/LxyErr",200,-100.0,100.0);
  AddOutput(hLxyOverLxyErr);
  
  hLz = new TH1F("Conversion Lz","Conversion Lz",200,-50.0,50.0);
  AddOutput(hLz);
  
  hLzOverLzErr = new TH1F("Conversion Lz/LzErr","Conversion Lz/LzErr",200,-100.0,100.0);
  AddOutput(hLzOverLzErr);
  
  hDxy = new TH1F("Conversion Dxy","Conversion Dxy",200,-5.0,5.0);
  AddOutput(hDxy);
  
  hSimMatchedLxy = new TH1F("SimMatched Conversion Lxy","SimMatched Conversion Lxy",100,-30.0,30.0);
  AddOutput(hSimMatchedLxy);
  
  hSimMatchedLz = new TH1F("SimMatched Conversion Lz","SimMatched Conversion Lz",200,-50.0,50.0);
  AddOutput(hSimMatchedLz);
  
  hSimMatchedDxy = new TH1F("SimMatched Conversion Dxy","SimMatched Conversion Dxy",200,-5.0,5.0);
  AddOutput(hSimMatchedDxy);

  hNTracks = new TH1F("Number of Tracks in the event", "Number of tracks in the event", 401, -0.5,400.5);
  AddOutput(hNTracks);
  
  hSimMatchedNTracks = new TH1F("SimMatched Number of Tracks in the event", "SimMatched Number of tracks in the event", 401, -0.5,400.5);
  AddOutput(hSimMatchedNTracks);
  
  hIsolation = new TH1F("Track Isolation of Conversion", "Track Isolation of Conversion", 201, -0.5,200.5);
  AddOutput(hIsolation);
  
  hSimMatchedIsolation = new TH1F("SimMatched Track Isolation of Conversion", "SimMatched Track Isolation of Conversion", 201, -0.5,200.5);
  AddOutput(hSimMatchedIsolation);

  hNConversions = new TH1F("Number of Conversions", "Number of Conversions",10,-0.5,10.5);
  AddOutput(hNConversions);

  hDoubleConvMass = new TH1F("Double Conversion Mass", "Double Conversion Mass", 1000, 0.0, 10.0);
  AddOutput(hDoubleConvMass);
  
  hConvGammaMass = new TH1F("Conversion plus Photon Mass", "Conversion plus Photon Mass", 1000, 0.0, 10.0);
  AddOutput(hConvGammaMass);
  
  hPiPiMass = new TH1F("Conversion two pion Mass", "Conversion two pion Mass", 1000,0.0,10.0);
  AddOutput(hPiPiMass);

  hAllConversionElectronPt = new TH1F("hAllConversionElectronPt","Electron Pt (GeV)",2000,0,200.0);
  AddOutput(hAllConversionElectronPt);
  
  hAllConversionElectronEta = new TH1F("hAllConversionElectronEta","Electron Eta",200,-5.0,5.0);
  AddOutput(hAllConversionElectronEta);
  
  hRemovedElectronPt = new TH1F("hRemovedElectronPt","Electron Pt (GeV)",2000,0,200.0);
  AddOutput(hRemovedElectronPt);
  
  hRemovedElectronEta = new TH1F("hRemovedElectronEta","Electron Eta",200,-5.0,5.0);
  AddOutput(hRemovedElectronEta);
  
  hAllConversionElectronPtWrongCharge = new TH1F("hAllConversionElectronPtWrongCharge","Electron Pt (GeV)",2000,0,200.0);
  AddOutput(hAllConversionElectronPtWrongCharge);
  
  hAllConversionElectronEtaWrongCharge = new TH1F("hAllConversionElectronEtaWrongCharge","Electron Eta",200,-5.0,5.0);
  AddOutput(hAllConversionElectronEtaWrongCharge);
  
  hRemovedElectronPtWrongCharge = new TH1F("hRemovedElectronPtWrongCharge","Electron Pt (GeV)",2000,0,200.0);
  AddOutput(hRemovedElectronPtWrongCharge);
  
  hRemovedElectronEtaWrongCharge = new TH1F("hRemovedElectronEtaWrongCharge","Electron Eta",200,-5.0,5.0);
  AddOutput(hRemovedElectronEtaWrongCharge);
  
  hAllConversionElectronChargeAssign = new TH1F("hAllConversionElectronChargeAssign","Electron Charge Assignment",3.0,-1.5,1.5);
  AddOutput(hAllConversionElectronChargeAssign);
  
  hRemovedElectronChargeAssign = new TH1F("hRemovedElectronChargeAssign","Electron Charge Assignment",3.0,-1.5,1.5);
  AddOutput(hRemovedElectronChargeAssign);
  
  hTrackAlgo = new TH1F("hTrackAlgo","Conversion Track Algorithm",101,-0.5,100.5);
  AddOutput(hTrackAlgo);
  
  hCharge = new TH1F("hCharge","Candidate Charge",5.0,-2.5,2.5);
  AddOutput(hCharge);
  
  hPt1Pt2 = new TH2F("hPt1Pt2","Pt1 vs Pt2", 500,0,10.0,500,0.0,10.0);
  AddOutput(hPt1Pt2);
  
  hClusVtxDiff = new TH1F("hClusVtxDiff","hClusVtxDiff",201,-0.5,200.5);
  AddOutput(hClusVtxDiff);
  
  hClusVtxQual = new TH1F("hClusVtxQual","hClusVtxQual",200,0.0,20);
  AddOutput(hClusVtxQual);
  
  fTimer.Start(1e8);

}

//__________________________________________________________________________________________________
void ConversionRemoval::Process()
{
  // Process entries of the tree. For this module, we just load the branch and fill the histograms.
  // --> Why is this done on the basis of the name? should be the pointer!

//   BitMask64 testMask;
//   testMask.SetBit(Track::TIB1S);
//   testMask.SetBit(Track::TIB2S);
//   testMask.SetBit(Track::TID1S);
//   testMask.SetBit(Track::TID2S);
//   testMask.SetBit(Track::TID3S);
//   testMask.SetBit(Track::TOB1S);
//   testMask.SetBit(Track::TOB2S);
//   testMask.SetBit(Track::TEC1S);
//   testMask.SetBit(Track::TEC2S);
//   testMask.SetBit(Track::TEC3S);
//   testMask.SetBit(Track::TEC4S);
//   testMask.SetBit(Track::TEC5S);
//   testMask.SetBit(Track::TEC6S);
//   testMask.SetBit(Track::TEC7S);
//   testMask.SetBit(Track::TEC8S);
//   testMask.SetBit(Track::TEC9S);
// 
//   const Long64_t *maskVal = reinterpret_cast<const Long64_t*>(testMask.Bits());
//   
//   printf("HitMaskVal = %lli\n",*maskVal);
//   
//   Track testTrack;
//   if (testMask == testTrack.StereoLayers())
//     printf("stereo mask match\n");
//   else
//     printf("stereo mask mismatch\n");
//   

  LoadBranch(fTrackName);
  LoadBranch("PrimaryVertexes");
  LoadBranch("EvtSelData");

  fCleanElectrons = GetObjThisEvt<ElectronCol>(fCleanElectronsName);
  ObjArray<Electron> *allConversionElectrons = new ObjArray<Electron>;
  AddObjThisEvt(allConversionElectrons,fAllConversionElectronsName);
  
  //printf("first loop\n")
  if (fComputeEff && (!fRequireCleanElectron || fCleanElectrons->GetEntries())) {
    LoadBranch(fMCPartName);
    for (UInt_t i=0; i<fCleanElectrons->GetEntries(); ++i) {
      const Electron *e = fCleanElectrons->At(i);
      if ( !fMatchSimConv && fMatchWe && SimWeMatch(e) )
        allConversionElectrons->Add(e);
      else if (fMatchSimConv && !fMatchWe && SimConversionMatch(e))
        allConversionElectrons->Add(e);
      else if (!fMatchSimConv && !fMatchWe)
        allConversionElectrons->Add(e);
    }
    
    for (UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
      const MCParticle* p = fMCParticles->At(i);
      const BitMask16 simFailed = FailedSimCuts(p);

      BitMask16 simPtFailed = simFailed;
      simPtFailed.ClearBit(eSimPt);
      if (!simPtFailed.NBitsSet())
        hSimPt->Fill(p->Pt());

      BitMask16 simEtaFailed = simFailed;
      simEtaFailed.ClearBit(eSimEta);
      if (!simEtaFailed.NBitsSet())
        hSimEta->Fill(p->Eta());
        
      BitMask16 simRhoFailed = simFailed;
      simRhoFailed.ClearBit(eSimRho);
      if ( !simRhoFailed.NBitsSet() ) {
          hSimConvRadius->Fill(p->DecayVertex().Rho());
          hSimConvPosition->Fill(p->DecayVertex().X(),p->DecayVertex().Y());
      }
      Double_t minPt=999;
      for (UInt_t j=0; j<p->NDaughters(); ++j) {
        const MCParticle *d = p->Daughter(j);
        if (TMath::Abs(d->PdgId())==11) {
          if (d->Pt()<minPt)
            minPt = d->Pt();
        }
      }
      
      BitMask16 simElectronPtFailed = simFailed;
      simElectronPtFailed.ClearBit(eSimElectronPt);
      if ( !simElectronPtFailed.NBitsSet() )
        hSimMinElectronPt->Fill(minPt);
        
    }
  }
  
  //printf("second loop\n")
  ObjArray<DecayParticle> *fakingConversions = new ObjArray<DecayParticle>;
  AddObjThisEvt(fakingConversions,fFakingConversionsName);
  if (!fRequireCleanElectron || allConversionElectrons->GetEntries()) {
    LoadBranch(fMvfConvName);
    for (UInt_t i=0; i<fMvfConversions->GetEntries(); ++i) {
      const DecayParticle* d = fMvfConversions->At(i);
      if (!fRequireCleanElectron || ElectronMatch(d,allConversionElectrons))
        fakingConversions->Add(d);
    }
  }
  
  
  //printf("third loop\n")
  ObjArray<DecayParticle> goodConversions;
  ObjArray<Electron> *removedElectrons = new ObjArray<Electron>;
  AddObjThisEvt(removedElectrons,fRemovedElectronsName);
  UInt_t nConversions = 0;
  for (UInt_t i=0; i<fakingConversions->GetEntries(); ++i) {
    //if (nConversions>=1) continue;
    const DecayParticle* c = fakingConversions->At(i);
    const BitMask16 failed = FailedCuts(c);
    
    // printf("Looking for sim match\n");
    const MCParticle *simPhoton = SimMatch(c);
    UInt_t nSimFailed=0;
    BitMask16 simFailed;
    //printf("Checking sim cuts\n");
    if (simPhoton)
      simFailed = FailedSimCuts(simPhoton);   
    
    //make N-2 plots for decay length
    BitMask16 decayLengthFailed = failed;
    decayLengthFailed.ClearBit(eRho);
    decayLengthFailed.ClearBit(eL);
    if (!decayLengthFailed.NBitsSet()) {
      hLxy->Fill(c->Lxy());
      hLxyOverLxyErr->Fill(c->Lxy()/c->LxyError());
      hLzOverLzErr->Fill(c->Lz()/c->LzError());      
      hLz->Fill(c->Lz());
      if (simPhoton && !simFailed.NBitsSet()) {
        hSimMatchedLxy->Fill(c->Lxy());
        hSimMatchedLz->Fill(c->Lz());
      }
    }
    
    //make N-1 plots
    //printf("Make N-1 plots\n");

    BitMask16 ptFailed = failed;
    ptFailed.ClearBit(ePt);
    if (!ptFailed.NBitsSet()) {
      hGammaPt->Fill(c->Pt());
      if (simPhoton) {
        BitMask16 simPtFailed = simFailed;
        simPtFailed.ClearBit(eSimPt);
        if (!simPtFailed.NBitsSet()) {
          hSimMatchedGammaPt->Fill(c->Pt());
          hSimMatchedSimGammaPt->Fill(simPhoton->Pt()); 
        }
      }
    }
    
    BitMask16 etaFailed = failed;
    etaFailed.ClearBit(eEta);
    if (!etaFailed.NBitsSet()) {
      hGammaEta->Fill(c->Eta());
      if (simPhoton) {
        BitMask16 simEtaFailed = simFailed;
        simEtaFailed.ClearBit(eSimEta);
        if (!simEtaFailed.NBitsSet()) {
          hSimMatchedGammaEta->Fill(c->Eta());
          hSimMatchedSimGammaEta->Fill(simPhoton->Eta());
        }
      }
    }
    
    BitMask16 rhoFailed = failed;
    rhoFailed.ClearBit(eRho);
    if ( !rhoFailed.NBitsSet() ) {
      hConversionRadius->Fill(c->Position().Rho());
      hConversionRPhi->Fill(c->Position().X(), c->Position().Y());
      hConversionRZ->Fill(c->Position().Z(), c->Position().Rho()*c->Position().Y()/TMath::Abs(c->Position().Y()));
      hConversionFoldedRZ->Fill(c->Position().Z(), c->Position().Rho());
      if (c->AbsEta()<1.2)
        hConversionRadiusBarrel->Fill(c->Position().Rho());
      if (simPhoton) {
        BitMask16 simRhoFailed = simFailed;
        simRhoFailed.ClearBit(eSimRho);
        if ( !simRhoFailed.NBitsSet() ) {
          hSimMatchedConvRadius->Fill(c->Position().Rho());
          hSimMatchedSimConvRadius->Fill(simPhoton->DecayVertex().Rho());
        }
      }    
    }
    
    BitMask16 dxyFailed = failed;
    dxyFailed.ClearBit(eDxy);
    if ( !dxyFailed.NBitsSet() ) {
      hDxy->Fill(c->Dxy());
      if (simPhoton && !simFailed.NBitsSet() )
        hSimMatchedDxy->Fill(c->Dxy());
    }
    
    BitMask16 probFailed = failed;
    probFailed.ClearBit(eProb);
    if (!probFailed.NBitsSet()) {
      hConvProb->Fill(c->Prob());
      hConvChi2->Fill(c->Chi2());
        if (simPhoton && !simFailed.NBitsSet() )
          hSimMatchedConvProb->Fill(c->Prob());
    }
    
    BitMask16 clusVtxFailed = failed;
    clusVtxFailed.ClearBit(eClusVtx);
    if (!clusVtxFailed.NBitsSet()) {
      hClusVtxDiff->Fill(fEvtSel->ClusVtxDiff());
      hClusVtxQual->Fill(fEvtSel->ClusVtxQual());
    }
    
    BitMask16 sharedHitsFailed = failed;
    sharedHitsFailed.ClearBit(eSharedHits);
    if (!sharedHitsFailed.NBitsSet()) {
      hNSharedHits->Fill(c->NSharedHits());
    }
    
    
    //printf("Loop over daughters\n");
    
    BitMask16 daughtersFailed = failed;
    daughtersFailed.ClearBit(eDaughters);
    if ( !daughtersFailed.NBitsSet() ) {
      Double_t minElectronPt = c->Daughter(0)->Pt();
      Double_t maxElectronPt = c->Daughter(1)->Pt();
      if (c->Daughter(1)->Pt() < c->Daughter(0)->Pt()) {
        minElectronPt = c->Daughter(1)->Pt();
        maxElectronPt = c->Daughter(0)->Pt();
      }
      
      const Particle *minPtElectron = c->Daughter(0);
      BitMask16 minPtElectronFailed;
      for (UInt_t j=0; j<1; ++j) {
        const DaughterData* d = c->DaughterDat(j);
        const StableData *s = dynamic_cast<const StableData*>(d);
        if (!s)
          Fatal("FailedElectronCuts", "DaughterData not a StableData");
        
        const ChargedParticle* dc = dynamic_cast<const ChargedParticle*>(c->Daughter(j));
        if (!dc)
          Fatal("FailedElectronCuts", "Daughter not a ChargedParticle");
        
        const Track* t = dc->Trk();
          
        const BitMask16 electronFailed = FailedElectronCuts(d);
        
        if (d->Pt() < minElectronPt) {
          minPtElectron = dc;
          minPtElectronFailed = electronFailed;
        }
        
        BitMask16 failedEPt = electronFailed;
        failedEPt.ClearBit(eEPt);
        if (!failedEPt.NBitsSet())
          hElectronPt->Fill(d->Pt());
          
        BitMask16 failedEEta = electronFailed;
        failedEEta.ClearBit(eEEta);
        if (!failedEEta.NBitsSet())
          hElectronEta->Fill(d->Eta());
          
        BitMask16 failedENHits = electronFailed;
        failedENHits.ClearBit(eENHits);
        if (!failedENHits.NBitsSet()) {
          hTrackNHits->Fill(t->NHits());
          if (1 )
            hSimMatchedTrackNHits->Fill(t->NHits());
        }
        
        BitMask16 failedENMissingHits = electronFailed;
        failedENMissingHits.ClearBit(eENMissingHits);
        if (!failedENMissingHits.NBitsSet()) {
          hTrackNMissingHits->Fill(s->NMissedHits());
        }
         
        BitMask16 failedENWrongHits = electronFailed;
        failedENWrongHits.ClearBit(eENWrongHits);
        if (!failedENWrongHits.NBitsSet()) {
          hTrackNWrongHits->Fill(s->NWrongHits());
        }
        
        BitMask16 failedEProb = electronFailed;
        failedEProb.ClearBit(eEProb);
        if (!failedEProb.NBitsSet()) {
          hTrackProb->Fill(t->Prob());
          if (1 ) {
            hSimMatchedTrackChi2->Fill(t->Chi2()/((Double_t)t->Ndof()));
            hSimMatchedTrackProb->Fill(t->Prob());
          }
        }
          
        BitMask16 failedEProbNHits = electronFailed;
        failedEProbNHits.ClearBit(eEProb);
        failedEProbNHits.ClearBit(eENHits);
        if ( !failedEProbNHits.NBitsSet() ) {
          hTrackNHitsProb->Fill(t->NHits(),t->RChi2());
          hTrackChi2->Fill(t->Chi2()/((Double_t)t->Ndof()));
          if (1)
            hSimMatchedTrackNHitsProb->Fill(t->NHits(),t->RChi2());
        }
        
        const mithep::Vertex *pv = fPrimaryVertexes->At(0);
        
        BitMask16 failedEDzPv = electronFailed;
        failedEDzPv.ClearBit(eDz);
        if (!failedEDzPv.NBitsSet()) {
          hTrackDzPv->Fill(t->Z0()-pv->Z());
        }
        
        if ( !electronFailed.NBitsSet()) {
          hTrackD0->Fill(t->D0());
        }
      }
      
      minPtElectronFailed.ClearBit(eEPt);
      if (!minPtElectronFailed.NBitsSet()) {
        hMinElectronPt->Fill(minElectronPt);
        hPt1Pt2->Fill(maxElectronPt,minElectronPt);
        if (simPhoton) {
          BitMask16 simElectronPtFailed = simFailed;
          simElectronPtFailed.ClearBit(eSimElectronPt);
          if ( !simElectronPtFailed.NBitsSet() ) {
//             if (!minPtElectron)
//               printf("no minptelectron");
            const ChargedParticle* minElectronC = dynamic_cast<const ChargedParticle*>(minPtElectron);
//             if (!minElectronC)
//               printf("no minelectronC\n");
            const Track *t = minElectronC->Trk();
//             if (!t)
//               printf("no t\n");
            if (t->MCPart())
              hSimMatchedMinElectronSimPt->Fill(t->MCPart()->Pt());
            hSimMatchedMinElectronPt->Fill(minElectronPt);
          }
        }
      }
    }
    
    //printf("Good conversions only\n");
    if ( !failed.NBitsSet() ) {
      goodConversions.Add(c);
      ++nConversions;
      //fill photon histograms
      hConversionZ->Fill(c->Position().Z());
      hConvDCotTheta->Fill(DCotTheta(c));
      hGammaMass->Fill(c->Mass());
      hCharge->Fill(c->Charge());
      hGammaPhi->Fill(c->Phi());
      
      for (UInt_t j=0; j<1; ++j) {
        const Particle* d = c->Daughter(j);
        const ChargedParticle* dc = dynamic_cast<const ChargedParticle*>(d);
        if (dc) {
          const Track* t = dc->Trk();
          hTrackAlgo->Fill(t->Algo());
        }
      }
      
      const Electron *removedElectron = ElectronMatch(c,allConversionElectrons);
      if (removedElectron) {
        removedElectrons->Add(removedElectron);
      }
      
      if (simPhoton && !simFailed.NBitsSet()) {

        //if (c->Position().Rho() < fRhoMax)
        hSimMatchedGammaMass->Fill(c->Mass());
        hSimMatchedConvResolution->Fill(c->Position().Rho() - simPhoton->DecayVertex().Rho());   
        hSimMatchedConvPhiRes->Fill(c->Position().Phi() - simPhoton->DecayVertex().Phi());

        for (UInt_t j=0; j<1; ++j) {
          const Particle* d = c->Daughter(j);
          const ChargedParticle* dc = dynamic_cast<const ChargedParticle*>(d);
          if (dc) {
            const Track* t = dc->Trk();
            hSimMatchedTrackD0->Fill(t->D0());
          }
        }
      }  
    }
  }
  hNConversions->Fill(nConversions);
  if (nConversions>0)
    printf("bx number = %i, nConversions = %i\n", GetEventHeader()->BunchCrossing(), nConversions);
  
  //loop over good conversions only
  for (UInt_t i=0; i<goodConversions.GetEntries(); ++i) {
    for (UInt_t j=i+1; j<goodConversions.GetEntries(); ++j) {
      CompositeParticle *p0 = new CompositeParticle;
      p0->AddDaughter(goodConversions.At(i));
      p0->AddDaughter(goodConversions.At(j));
      hPiPiMass->Fill(p0->Mass());
      delete p0;
    }
  }
  
  
  for (UInt_t i=0; i<allConversionElectrons->GetEntries(); ++i) {
    const Electron *e = allConversionElectrons->At(i);
    hAllConversionElectronPt->Fill(e->Pt());
    hAllConversionElectronEta->Fill(e->Eta());
    const MCParticle *mcElectron = SimWeMatch(e);
    if (mcElectron) {
      const Double_t chargeAssign = e->Charge()*mcElectron->Charge();
      hAllConversionElectronChargeAssign->Fill(chargeAssign);
      if (chargeAssign==-1.0) {
        hAllConversionElectronPtWrongCharge->Fill(e->Pt());
        hAllConversionElectronEtaWrongCharge->Fill(e->Eta());
      }
    }
    else
      hAllConversionElectronChargeAssign->Fill(0.0);
  }
  
  for (UInt_t i=0; i<removedElectrons->GetEntries(); ++i) {
    const Electron *e = removedElectrons->At(i);
    hRemovedElectronPt->Fill(e->Pt());
    hRemovedElectronEta->Fill(e->Eta());
    const MCParticle *mcElectron = SimWeMatch(e);
    if (mcElectron) {
      const Double_t chargeAssign = e->Charge()*mcElectron->Charge();
      hRemovedElectronChargeAssign->Fill(chargeAssign);
      if (chargeAssign==-1.0) {
        hRemovedElectronPtWrongCharge->Fill(e->Pt());
        hRemovedElectronEtaWrongCharge->Fill(e->Eta());
      }
    }
    else
      hRemovedElectronChargeAssign->Fill(0.0);
  }


}

//__________________________________________________________________________________________________
void ConversionRemoval::SlaveTerminate()
{
  // Run finishing code on the computer (slave) that did the analysis. For this module, we dont do
  // anything here.
}

//__________________________________________________________________________________________________
void ConversionRemoval::Terminate()
{

  fTimer.Stop();
  printf("Total Analysis Time = %f\n",fTimer.CpuTime());


  TCanvas* c1 = new TCanvas();
  hConversionRadius->Draw();
  
  TCanvas* c2 = new TCanvas();
  hConversionRPhi->Draw("col");
  //c2->SetLogz();
  
  // Run finishing code on the client computer. For this module, we dont do anything here.
}

BitMask16 ConversionRemoval::FailedCuts(const DecayParticle* c) {


        //if (p->GetVertex()->Prob() < 0.005)
        //      return false;
        UInt_t nFailed = 0;
        
        BitMask16 failed;
        
        if ( c->NDaughters() != 2)
          failed.SetBit(eNDaughters);
          
        if ( c->Pt() < fPtMin )
          failed.SetBit(ePt);
          
        if ( !PassRho(c->Position().Rho()) )
          failed.SetBit(eRho);
          
        if ( TMath::Abs(c->Eta()) > fEtaMax )
          failed.SetBit(eEta); 
          
        if ( c->Lxy() < fLxyMin || c->Lz() < fLzMin )
          failed.SetBit(eL);
                    
        if ( TMath::Abs(c->Dxy()) > fAbsDxyMax )
          failed.SetBit(eDxy);
                 
        if ( TMath::Abs(c->Charge()) != 0 )
          failed.SetBit(eCharge);
        
//         if (GetEventHeader()->BunchCrossing()!=51 && GetEventHeader()->BunchCrossing()!=2724)
//           failed.SetBit(eBx);
        
//         if (fTracks->GetEntries()>150)
//           failed.SetBit(eBx);
        
//         Double_t absPhi = TMath::Abs(c->Phi());
//         if (!( absPhi<(TMath::Pi()/4.0) || absPhi>(3.0*TMath::Pi()/4.0)))
//           failed.SetBit(eBx);
//       
        if (c->NSharedHits() > fMaxNSharedHits || c->NSharedHits() < fMinNSharedHits)
          failed.SetBit(eSharedHits);

        const mithep::Vertex *pv = fPrimaryVertexes->At(0);
        if (fRequireGoodPv && (pv->NTracksFit()<2 || pv->Prob()<1e-6 || pv->Z()<-12.5 || pv->Z()>12.5) )
          failed.SetBit(eBx);
      
        if ( fEvtSel->ClusVtxQual() < fMinClusterVtxQual )
          failed.SetBit(eClusVtx);
        
        Bool_t failedElectrons=false;
        for (UInt_t i=0; i<c->NDaughters(); ++i) {
          const DaughterData* d = c->DaughterDat(i);
          const Track *trk = dynamic_cast<const ChargedParticle*>
                       (c->Daughter(i))->Trk();
          Bool_t trackQuality = ( !trk->IsGsf() || !fTrackQualityFirstOnly );
          BitMask16 failedElectronCuts = FailedElectronCuts(d, trackQuality);
          if (failedElectronCuts.NBitsSet())
            failedElectrons=true;
        }
        if (failedElectrons)
          failed.SetBit(eDaughters);
        
        if ( c->Prob() < fProbMin )
          failed.SetBit(eProb);

        return failed;
//         LoadBranch(fMvfConvUnconstrainedName);
//         const DecayParticle *unconstrained = MatchingConversion(c, fMvfConversionsUnconstrained);
// 
//         if (!unconstrained)
//           nFailed++;
// 
//         if ( TwoPionMass(unconstrained) > 0.4 )
//           nFailed++;
        //hConvDCotTheta->Fill(dCotTheta);

}

Bool_t ConversionRemoval::PassConvCuts(const Conversion* p) {
        if (p->NDaughters() != 2)
        return false;

        Double_t convRadius = p->DecayVertex().Position().Rho();
        if (convRadius < 2.0)
        //if (convRadius < 2.0 || convRadius > 20.0)
        //if (convRadius < 20.0)
                return false;
        
        const Electron* pElectron1 = (Electron*)p->Daughter(0);
        const Electron* pElectron2 = (Electron*)p->Daughter(1);
        
        
        //Double_t deltaPhi = TMath::ACos((pElectron1->Px()*pElectron2->Px() + pElectron1->Py()*pElectron2->Py() + pElectron1->Pz()*pElectron2->Pz())/(pElectron1->Mom().P()*pElectron2->Mom().P()));
        
        //if (deltaPhi>0.5)
        //      return false;
       // if (fabs(p->DCotTheta())>0.2)
       if (fabs(p->DCotTheta())>0.75)
          return false;
                
        if (p->Charge()!=0)
                return false;
                
                
/*        for (Int_t i=0; i<p->NDaughters(); i++) {
          const Electron* e = p->Daughter(i);
          if (!PassElectronCuts(e))
            return false;
        }  */ 
        //if (p->Mass()<0.6)
        //      return false;

        //if (p->GetVertex().Prob() < 0.2)// || p->GetVertex().Prob() > 0.9)
        //      return false;
        
        return true;
}

BitMask16 ConversionRemoval::FailedElectronCuts(const DaughterData* d, Bool_t trackQuality) {
        //if (p->Pt() < 5.0)
        //  return false;
        
        BitMask16 failed;
        
        if (d->Pt() < fElectronPtMin)
          failed.SetBit(eEPt);
          
        if (TMath::Abs(d->Eta()) > fElectronEtaMax)
          failed.SetBit(eEEta);
        
        const StableData *s = dynamic_cast<const StableData*>(d);

        if (!s)
          Fatal("FailedElectronCuts", "DaughterData not a StableData");

        const ChargedParticle* c = dynamic_cast<const ChargedParticle*>(s->Original());
        
        if (!c)
          Fatal("FailedElectronCuts", "Daughter not a ChargedParticle");
          
        const Track* t = c->Trk();
          
        if (!t)
          Fatal("FailedElectronCuts", "Daughter has no track");
        
        if (s->NMissedHits() > fMissedHitsMax)
          failed.SetBit(eENMissingHits);

        if (s->NWrongHits() > fWrongHitsMax)
          failed.SetBit(eENWrongHits);
        
        if (trackQuality && t->NHits()<fNHitsMin)
          failed.SetBit(eENHits);
  
        if (trackQuality && t->NStereoHits()<fNStereoHitsMin)
          failed.SetBit(eENStereoHits);
        
        if (trackQuality && t->NPixelHits()<fNPixelHitsMin)
          failed.SetBit(eENPixelHits);
        
        if (trackQuality && t->Prob()<fTrackProbMin)
          failed.SetBit(eEProb);
          
        //if (trackQuality && fExcludePXB1 && t->Hit(Track::PXB1))
        if (trackQuality && fExcludePXB1 && (t->Hit(Track::PXB1)||t->Hit(Track::PXF1)))
          failed.SetBit(eEPxb1);

        const mithep::Vertex *pv = fPrimaryVertexes->At(0);
        
        if (trackQuality && (TMath::Abs(t->Z0()-pv->Z()) > fMaxTrackDzPv) )
          failed.SetBit(eDz);
        
        return failed;
}

Bool_t ConversionRemoval::PassTrackCuts(const Track* t) {

 // if ( t->Prob()<fTrackProbMin && t->NHits()<fNHitsMin )

//   if (t->NStereoHits() < 2)
//     return false;
    
  return false;

}

BitMask16 ConversionRemoval::FailedSimCuts(const MCParticle* p) {
        
        BitMask16 failed;
        
        if (p->PdgId()!=22)
                failed.SetBit(eSimPdg);
        if (TMath::Abs(p->Eta())>fEtaMax)
                failed.SetBit(eSimEta);
        //if (p->Pt()<5.0 || p->Pt()>25.0)
        //if (p->Pt()<25.0)
        if (p->Pt()<fPtMin)
                failed.SetBit(eSimPt);
                
        Double_t simRad = p->DecayVertex().Rho();
        //if (simRad<2.0 || simRad>128.0)
        //if (simRad<2.0 || simRad>20.0)
        //if (simRad<20.0 || simRad>128.0)
        if ( !PassRho(simRad) )
          failed.SetBit(eSimRho);
          
        Int_t nEPlus = 0;
        Int_t nEMinus = 0;
        Bool_t failedElectronPt=false;
        Bool_t failedElectronEta=false;
        for (UInt_t j=0; j<p->NDaughters(); j++) {
                const MCParticle* pDaughter = p->Daughter(j);
                if (pDaughter->PdgId()==11)
                        nEPlus++;
                if (pDaughter->PdgId()==-11)
                        nEMinus++;
                if (TMath::Abs(pDaughter->PdgId())==11) {
                  if (pDaughter->Pt() < fElectronPtMin)
                    failedElectronPt = true;
                  if (TMath::Abs(pDaughter->Eta()) > fElectronEtaMax)
                    failedElectronEta=true;
                }
        }
        if (failedElectronPt)
          failed.SetBit(eSimElectronPt);
        if (failedElectronEta)
          failed.SetBit(eSimElectronEta);
        
        if (nEPlus!=1 || nEMinus!=1)
          failed.SetBit(eSimDaughters);
        
        //get rid of all conversions which either don't produce a clean electron, or which
        //were produced from a real primary electron
        //if (fRequireCleanElectron && ( !ElectronMatch(p) || GenElectronMatch(p) ) )
        if (fRequireCleanElectron && !ElectronMatch(p))
          failed.SetBit(eCleanElectron);
          
        return failed;
}

const Electron* ConversionRemoval::ElectronMatch(const MCParticle *p) {
//   if (p->HasMother(22))
//     return 0;
  const Electron *eMatch = 0;
  if (fHitBasedMatching) {
    for (UInt_t i=0; i<fCleanElectrons->GetEntries(); ++i) {
      const Electron *e = fCleanElectrons->At(i);
      const MCParticle *ep = e->Trk()->MCPart();
      if ( ep && ep->DistinctMother()==p)
        return e;
    }
  }
  else {
    Double_t minDeltaR = 0.3;
    for (UInt_t i=0; i<fCleanElectrons->GetEntries(); ++i) {
      const Electron *e = fCleanElectrons->At(i);
      Double_t deltaR = MathUtils::DeltaR(*e,*p);
      if (deltaR < minDeltaR) {
        minDeltaR = deltaR;
        eMatch = e;
      }
    }
  }
  
  return eMatch;

}

const MCParticle* ConversionRemoval::GenElectronMatch(const MCParticle *p) {
  const MCParticle *mother = p;
  while (p) {
    if ( p->AbsPdgId()==11 && p->IsGenerated() )
      return p;
    else
      p = p->Mother();
  }
  return 0;
}

const MCParticle* ConversionRemoval::SimMatch(const DecayParticle* p) {
  if ( p->NDaughters()!=2 )
    return 0;
    
  const ChargedParticle *electron1 = dynamic_cast<const ChargedParticle*>(p->Daughter(0));
  const ChargedParticle *electron2 = dynamic_cast<const ChargedParticle*>(p->Daughter(1));
  
  if (!electron1 || !electron2)
    return 0;
    
  const MCParticle *simElectron1 = 0;
  const MCParticle *simElectron2 = 0;    
  if (fHitBasedMatching) {
    simElectron1 = electron1->Trk()->MCPart();
    simElectron2 = electron2->Trk()->MCPart();
  }
  else {
    LoadBranch(fMCPartName);
    for (UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
      const MCParticle *mcp = fMCParticles->At(i);
      if (mcp->PdgId()==22 && mcp->DecayVertex().Rho()<80.0 && MathUtils::DeltaR(*mcp,*p) < 0.3)
        return mcp;
    }
    return 0;
  }
    
  if (simElectron1) {
    hTrackSimMatchType->Fill(simElectron1->PdgId());
    if (simElectron1->HasMother())
      hParentSimMatchType->Fill(simElectron1->Mother()->PdgId());
  }
  
  if (simElectron2) {
    hTrackSimMatchType->Fill(simElectron2->PdgId());
    if (simElectron2->HasMother())
      hParentSimMatchType->Fill(simElectron2->Mother()->PdgId());
  }
    
  
  if (!simElectron1) {
    const Track *t = electron1->Trk();
    hUnMatchedTrackChi2->Fill(t->Chi2()/((Double_t)t->Ndof()));
    hUnMatchedTrackProb->Fill(t->Prob());
    hUnMatchedTrackNHits->Fill(t->NHits());
    hUnMatchedTrackNHitsProb->Fill(t->NHits(),t->RChi2());
    hUnMatchedTrackD0->Fill(t->D0());
  }
  
  if (!simElectron2) {
    const Track *t = electron2->Trk();
    hUnMatchedTrackChi2->Fill(t->Chi2()/((Double_t)t->Ndof()));
    hUnMatchedTrackProb->Fill(t->Prob());
    hUnMatchedTrackNHits->Fill(t->NHits());
    hUnMatchedTrackNHitsProb->Fill(t->NHits(),t->RChi2());
    hUnMatchedTrackD0->Fill(t->D0());
  }
  
  if (!simElectron1 || !simElectron2) {
    //printf("unmatched electron\n");
    return 0;
  }
    
  
  if ( TMath::Abs(simElectron1->PdgId())!=11  || TMath::Abs(simElectron2->PdgId())!=11 ) {
    //printf("not electrons\n");
    return 0;
  }
      
  const MCParticle *simPhoton = simElectron1->DistinctMother();
  
  if (!simPhoton) {
   // printf("no photon\n");
  //  hParentSimMatchType->Fill(0);
    return 0;
  }
    
  if (simElectron2->DistinctMother()!=simPhoton) {
    //printf("different mothers\n");
   // hParentSimMatchType->Fill(0);
    return 0;
  }
    
  if (simPhoton->PdgId()!=22) {
    //printf("not a photon\n");
    return 0; 
  }

 // printf("found common mother with pdg=%i, daughter pdgs: %i, %i\n", simPhoton->PdgId(), simElectron1->PdgId(),simElectron2->PdgId());
  
 // hParentSimMatchType->Fill(simPhoton->PdgId());
    
  return simPhoton;
  

}

const MCParticle* ConversionRemoval::SimMatch(const Track* t) {
    
  const MCParticle *simElectron1 = 0;
  //const MCParticle *simElectron2 = 0;    
  if (fHitBasedMatching) {
    simElectron1 = t->MCPart();
  }
  else {
    LoadBranch(fMCPartName);
    //Double_t minDeltaR = 999.9;
    for (UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
      const MCParticle *p = fMCParticles->At(i);
      if (p->PdgId()==22 && p->DecayVertex().Rho()<80.0 && MathUtils::DeltaR(*p,*t) < 0.3)
        return p;
    }
  }
  
  if (!simElectron1) {
    //printf("unmatched electron\n");
    return 0;
  }
    
  
  if ( TMath::Abs(simElectron1->PdgId())!=11 ) {
    //printf("not electrons\n");
    return 0;
  }
      
  const MCParticle *simPhoton = simElectron1->DistinctMother();
  
  if (!simPhoton) {
   // printf("no photon\n");
  //  hParentSimMatchType->Fill(0);
    return 0;
  }
    
  if (simPhoton->PdgId()!=22) {
    //printf("not a photon\n");
    return 0; 
  }

 // printf("found common mother with pdg=%i, daughter pdgs: %i, %i\n", simPhoton->PdgId(), simElectron1->PdgId(),simElectron2->PdgId());
  
 // hParentSimMatchType->Fill(simPhoton->PdgId());
    
  return simPhoton;
  

}

UInt_t ConversionRemoval::NTracks(const TrackCol *col)
{
  UInt_t ntracks=0;
  for (UInt_t i=0; i<col->GetEntries(); ++i) {
    const Track *t = col->At(i);
    if ( PassTrackCuts(t) )
      ntracks++;
  }
  
  return ntracks;
}

// Double_t ConversionRemoval::TrackPtIsolation(const DecayParticle *p, Double_t r)
// {
//   TObjArray conversionTracks;
//   for (UInt_t i=0; i<p->NDaughters(); ++i) {
//     const ChargedParticle* d = dynamic_cast<const ChargedParticle*>(p->Daughter(i));
//     if (d) {
//       const Track *t = d->Trk();
//       if (t)
//         conversionTracks.Add((TObject*)t);
//     }
//   }
//   
//   
//   LoadBranch(fTrackName);
//   Double_t sumPt=0.0;
//   for (UInt_t i=0; i<fTracks->GetEntries(); ++i) {
//     const Track *t = fTracks->At(i);
//     if (conversionTracks.IndexOf(t)==-1) {
//       Double_t deltaR = MathUtils::DeltaR(p->Phi(),p->Eta(),t->Phi(),t->Eta());
//       if (deltaR<r)
//         sumPt = sumPt + 1.0;//sumPt += t->Pt();
//     }
//   }
//    
//   return sumPt; 
// 
// }

Bool_t ConversionRemoval::PassPhotonCuts(const Photon *p)
{
  if (p->IsConverted())
    return false;
  
  return true;
}

Double_t ConversionRemoval::DCotTheta(const DecayParticle *c)
{
  if (c->NDaughters()!=2)
    return 0.0;

  Double_t theta1 =  c->DaughterDat(0)->Theta();
  Double_t theta2 =  c->DaughterDat(1)->Theta();
  
  Double_t dCotTheta = 1./TMath::Tan(theta1) - 1./TMath::Tan(theta2);
  //Double_t dCotTheta = c->DaughterDat(0)->Pt() - c->DaughterDat(1)->Pt();  

  if (c->DaughterDat(0)->Charge()==1)
    return dCotTheta;
  else
    return (-dCotTheta);
}

Double_t ConversionRemoval::DPhi(const DecayParticle *c)
{
  if (c->NDaughters()!=2)
    return 0.0;

  ThreeVectorC mom1 =  ((const StableData*)c->DaughterDat(0))->ThreeMom();
  ThreeVectorC mom2 =  ((const StableData*)c->DaughterDat(1))->ThreeMom();
  
  Double_t dPhi = TMath::ACos(mom1.Dot(mom2)/(mom1.R()*mom2.R()));
  
  return dPhi;
}

Double_t ConversionRemoval::TwoPionMass(const DecayParticle *c)
{
  if (c->NDaughters()!=2)
    return 0.0;

  StableParticle *pi1 = new StableParticle(211, (Track*)((ChargedParticle*)c->Daughter(0))->Trk());
  StableParticle *pi2 = new StableParticle(211, (Track*)((ChargedParticle*)c->Daughter(1))->Trk());
  
  StableData *pi1d = new StableData(pi1, ((StableData*)c->DaughterDat(0))->ThreeMom());
  StableData *pi2d = new StableData(pi2, ((StableData*)c->DaughterDat(1))->ThreeMom());
  
  CompositeParticle *meson = new CompositeParticle();
  meson->AddDaughter(pi1d);
  meson->AddDaughter(pi2d);
    
  Double_t mesonMass = meson->Mass();
  
  delete meson;
  delete pi1d;
  delete pi2d;
  delete pi1;
  delete pi2;
  
  return mesonMass;
}

const DecayParticle *ConversionRemoval::MatchingConversion(const DecayParticle *c, const DecayParticleCol *col) const
{
  for (UInt_t i=0; i<col->GetEntries(); ++i) {
    const DecayParticle *p = col->At(i);
    if ( p->HasSameDaughters(c) )
      return c;
  }

  return 0;
}

const Electron *ConversionRemoval::MatchingElectron(const ChargedParticle *c, const ElectronCol *col) const
{
   for (UInt_t i=0; i<col->GetEntries(); ++i) {
     const Electron *e = col->At(i);

    if ( c->Trk() == e->Trk() )
      return e;
   }
  
  return 0;
}

const Electron *ConversionRemoval::MatchingElectron(const Electron *ein, const ElectronCol *col) const
{
   for (UInt_t i=0; i<col->GetEntries(); ++i) {
     const Electron *e = col->At(i);

    if ( e->SCluster() == ein->SCluster() )
      return e;
   }
  
  return 0;
}

void ConversionRemoval::AddRhoRange(Double_t lower, Double_t upper)
{
  fRhoLbs.push_back(lower);
  fRhoUbs.push_back(upper);
}

Bool_t ConversionRemoval::PassRho(Double_t rho)
{
  for (UInt_t i=0; i<fRhoLbs.size(); ++i) {
    if ( rho>fRhoLbs.at(i) && rho<fRhoUbs.at(i) )
      return true;
  }
  
  return false;
}

const MCParticle *ConversionRemoval::SimConversionMatch(const Electron *e) {

  //LoadBranch(fMCPartName);

  if (fHitBasedMatching) {
    const MCParticle *p = e->Trk()->MCPart();
    if (p && p->DistinctMother() && !FailedSimCuts(p->DistinctMother()).NBitsSet())
      return p->DistinctMother();
  
  }
  else {
    for (UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
      const MCParticle *p = fMCParticles->At(i);
      //if (MathUtils::DeltaPhi(*e,*p) < 0.3 && !FailedSimCuts(p).NBitsSet())
      if ( MathUtils::DeltaPhi(*e,*p) < 0.3 && p->AbsPdgId()==11 && p->DistinctMother() && !FailedSimCuts(p->DistinctMother()).NBitsSet())
        return p->DistinctMother();
    }
  }
  
  return 0;

}


const MCParticle *ConversionRemoval::SimWeMatch(const Electron *e) {

  if (fHitBasedMatching) {
    const MCParticle *p = e->Trk()->MCPart();
    if (p && p->AbsPdgId()==11 && p->DistinctMother() && (p->DistinctMother()->AbsPdgId()==24 || p->DistinctMother()->AbsPdgId()==23) )
      return p;  
  }
  else {
    for (UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
      const MCParticle *p = fMCParticles->At(i);
      if (MathUtils::DeltaPhi(*e,*p) < 0.3 && p->AbsPdgId()==11 && p->DistinctMother() && (p->DistinctMother()->AbsPdgId()==24 || p->DistinctMother()->AbsPdgId()==23) )
        return p;
    }
  }
  
  return 0;
}

const Electron *ConversionRemoval::ElectronMatch(const DecayParticle *d, const ElectronCol *c) {

  if (d->NDaughters()<2)
    return 0;
    
  const ChargedParticle *cp = dynamic_cast<const ChargedParticle*>(d->Daughter(1));
  if (!cp)
    return 0;
    
  const Track *t = cp->Trk();
  if (!t)
    return 0;
    
  for (UInt_t i=0; i<c->GetEntries(); ++i) {
    const Electron *e = c->At(i);
    if (e->Trk()==t)
      return e;
  }
  
  return 0;
  
}

Revision:	1.1
Committed:	Wed Dec 16 18:10:38 2009 UTC (15 years, 4 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
Log Message:	updated MitConversions