Mods/src/MvfConversions.cc

// $Id: MvfConversions.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/MvfConversions.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::MvfConversions)

//__________________________________________________________________________________________________
MvfConversions::MvfConversions(const char *name, const char *title) :
  BaseMod    (name,title),
  fPtMin   (0),
  fEtaMax  (0),
  fProbMin (0),
  fExcludePXB1(false),
  fMassMax(0),
  fMissedHitsMax(0),
  fWrongHitsMax(0),
  fElectronPtMin(0),
  fElectronEtaMax(0),
  fComputeEff(kFALSE),
  fHitBasedMatching(kFALSE),
  fRequireCleanElectron(kFALSE),
  fTrackQualityFirstOnly(kFALSE),
  fTracks     (0),
  fGsfTracks  (0),
  fInOutTracks (0),
  fOutInTracks (0),
  fElectrons  (0),
  fTrkElectrons(0),
  fPhotons (0),
  fMCParticles (0),
  fMvfConversions (0),
  fMvfConversionsUnconstrained (0),
  fCleanElectrons(0),
  fMCPartName (Names::gkMCPartBrn),
  fTrackName  (Names::gkTrackBrn),
  fConvElectronName   ("Electrons"),
  fPhotonName   (Names::gkPhotonBrn),
  fMvfConvName("MvfConversions"),
  fMvfConvUnconstrainedName("MvfConversionsUnconstrained"),
  fElectronName(Names::gkElectronBrn),
  fTrkElectronName("TrackerElectrons"),
  fGoodElectronsName("ConversionElectrons"),
  fGoodTwoClusterElectronsName("TwoClusterConversionElectrons"),
  fBadElectronsName("ConversionRemovedElectrons"),
  hConversionRadius   (0),
  hConversionZ(0),
  hConversionRPhi (0),
  hConversionRZ(0),
  hGammaPt (0),
  hGammaEta (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),
  hTrackNHitsProb     (0),
  hTrackD0 (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) 
{
  // Constructor.
}

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

//__________________________________________________________________________________________________
void MvfConversions::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);

  // 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);
  
  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);
  
  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);
  
  hGammaMass = new TH1F("Photon Invariant Mass","Photon Invariant Mass",200,0,20.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)",400,0,200.0);
  AddOutput(hMinElectronPt);
  
  hSimMatchedMinElectronPt = new TH1F("Sim Matched Min Electron Pt","Sim Matched Min Electron Pt (GeV)",400,0,200.0);
  AddOutput(hSimMatchedMinElectronPt);
  
  hSimMatchedMinElectronSimPt = new TH1F("Sim Matched Min Electron Sim Pt","Sim Matched Min Electron Sim Pt (GeV)",400,0,200.0);
  AddOutput(hSimMatchedMinElectronSimPt);
  
  hSimMinElectronPt = new TH1F("Sim Min Electron Pt","Sim Min Electron Pt (GeV)",400,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",400,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",400,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",400,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)",400,0,200.0);
  AddOutput(hSimMatchedElectronPt);
  
  hSimMatchedGammaPt = new TH1F("Sim Matched Photon Pt","Sim Matched Photon Pt (GeV)",400,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)",400,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);
  
  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);
  
  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);

  fTimer.Start(1e8);

}

//__________________________________________________________________________________________________
void MvfConversions::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");
//   

  if (fRequireCleanElectron)
    fCleanElectrons = GetObjThisEvt<ElectronCol>(fCleanElectronsName);
  
  //printf("Loop over sim conversions\n");
  if (fComputeEff) {
    LoadBranch(fMCPartName);
    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);
        
      if (fRequireCleanElectron && !simFailed.NBitsSet()) {
        printf("Found good sim conversion, printing mother chain:\n");
        const MCParticle *pm = p;
        while (pm) {
          pm->Print();
          pm = pm->Mother();
        }
        const Electron *me = ElectronMatch(p);
        if (me) {
          printf("Matching Electron:\n");
          me->Print();
        }
        const MCParticle *mce = GenElectronMatch(p);
        if (mce) {
          printf("Matching Gen Electron:\n");
          mce->Print();
        }
      }
    
    }
  }
 
//   Track* trackTest = new Track();
//   Electron* electronTest = new Electron();
//   electronTest->SetTrackerTrk(trackTest);
//   
//   Track *trackCopy = new Track(*trackTest);
//   Electron* electronCopy = new Electron(*electronTest);
//   
//   const Track *electronTrack = electronTest->TrackerTrk();
//   const Track *electronCopyTrack = electronCopy->TrackerTrk();
//   
//   Electron *secondElectron = new Electron();
//   secondElectron->SetTrackerTrk(trackCopy);
//   
//   const Track *secondElectronTrack = secondElectron->TrackerTrk();  
//   
//   if (electronTrack==trackTest)
//     printf ("electronTrack==trackTest\n");
//     
//   if (electronTrack==trackCopy)
//     printf("electronTrack==trackCopy\n");
//     
//   if (electronCopyTrack==trackTest)
//     printf("electronCopyTrack==trackTest\n");
//   
//     
//   if (electronCopyTrack==trackCopy)
//     printf("electronCopyTrack==trackCopy\n");
//     
//   if (secondElectronTrack==trackTest)
//     printf("secondElectronTrack==trackTest\n");
//     
//   if (secondElectronTrack==trackCopy)
//     printf("secondElectronTrack==trackCopy\n");


  //printf("Loop over reco conversions\n");
  LoadBranch(fMvfConvName);
  ObjArray<DecayParticle> goodConversions;
  ObjArray<Electron> *goodElectrons = new ObjArray<Electron>;
  ObjArray<Electron> *goodTwoClusterElectrons = new ObjArray<Electron>;
  ObjArray<Electron> *badElectrons = new ObjArray<Electron>;
  UInt_t nConversions = 0;
  for (UInt_t i=0; i<fMvfConversions->GetEntries(); ++i) {
    const DecayParticle* c = fMvfConversions->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);        
  
    //printf("getting on with life\n");
    
    //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()));
      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());
    }
    
    //printf("Loop over daughters\n");
    
    BitMask16 daughtersFailed = failed;
    daughtersFailed.ClearBit(eDaughters);
    if ( !daughtersFailed.NBitsSet() ) {
      Double_t minElectronPt = c->Daughter(0)->Pt();
      const Particle *minPtElectron = c->Daughter(0);
      BitMask16 minPtElectronFailed;
      for (UInt_t j=0; j<c->NDaughters(); ++j) {
        const DaughterData* d = c->DaughterDat(j);
        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) {
          minElectronPt = d->Pt();
          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 (simPhoton && !simFailed.NBitsSet() )
            hSimMatchedTrackNHits->Fill(t->NHits());
        }
        
        BitMask16 failedEProb = electronFailed;
        failedEProb.ClearBit(eEProb);
        if (!failedEProb.NBitsSet()) {
          hTrackProb->Fill(t->Prob());
          if (simPhoton && !simFailed.NBitsSet() ) {
            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 (simPhoton && !simFailed.NBitsSet() )
            hSimMatchedTrackNHitsProb->Fill(t->NHits(),t->RChi2());
        }
        
        if ( !electronFailed.NBitsSet()) {
          hTrackD0->Fill(t->D0());
        }
      }
      
      minPtElectronFailed.ClearBit(eEPt);
      if (!minPtElectronFailed.NBitsSet()) {
        hMinElectronPt->Fill(minElectronPt);
        if (simPhoton) {
          BitMask16 simMinPtElectronFailed = simFailed;
          simMinPtElectronFailed.ClearBit(eSimElectronPt);
          if (!simMinPtElectronFailed.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());
      
      //LoadBranch(fTrkElectronName);
      fTrkElectrons = new ElectronOArr;
      const SuperCluster* cluster1 = static_cast<const ChargedParticle*>(c->Daughter(0))->Trk()->SCluster();
      const SuperCluster* cluster2 = static_cast<const ChargedParticle*>(c->Daughter(1))->Trk()->SCluster();
      for (UInt_t j=0; j<c->NDaughters(); ++j) {
        const ChargedParticle* dc = dynamic_cast<const ChargedParticle*>(c->Daughter(j));
        const Electron *e = MatchingElectron(dc,fTrkElectrons);
        if (e) {
          //printf("Matching Electron added\n");
          if (cluster1 && cluster1!=cluster2)
            if (!goodTwoClusterElectrons->HasObject(e))
              goodTwoClusterElectrons->Add(e);
          
          if ( !goodElectrons->HasObject(e) )
            goodElectrons->Add(e);
        }
      }
      
      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<c->NDaughters(); ++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);
  
  //Fill bad electron collection (all electrons, but with good electrons from conversions removed)
  //LoadBranch(fElectronName);
  fElectrons = new ElectronOArr;
  for (UInt_t i = 0; i<fElectrons->GetEntries(); ++i) {
    const Electron *e = fElectrons->At(i);
    const Electron *goodElectronMatch = MatchingElectron(e,goodElectrons);
    if (!goodElectronMatch)
      badElectrons->Add(e);
  }
  
  AddObjThisEvt(goodElectrons,  fGoodElectronsName.Data());
  AddObjThisEvt(goodTwoClusterElectrons,  fGoodTwoClusterElectronsName.Data());
  AddObjThisEvt(badElectrons,   fBadElectronsName.Data());
  
//   LoadBranch(fPhotonName);
  for (UInt_t i=0; i<goodConversions.GetEntries(); i++) {
    const DecayParticle *conv1 = goodConversions.At(i);
    for (UInt_t j=i+1; j<goodConversions.GetEntries(); j++) {
      DecayParticle *conv2 = goodConversions.At(j);
      if (!conv2->HasCommonDaughter(conv1)) {
        CompositeParticle pi0;
        pi0.AddDaughter(conv1);
        pi0.AddDaughter(conv2);
        hDoubleConvMass->Fill(pi0.Mass());
      }
    }
//     for (UInt_t j=0; j<fPhotons->GetEntries(); j++) {
//       Photon *photon = fPhotons->At(j);
//       if (PassPhotonCuts(photon)) {
//         CompositeParticle pi0;
//         pi0.AddDaughter(conv1);
//         pi0.AddDaughter(photon);
//         hConvGammaMass->Fill(pi0.Mass());
//       }
//     }
   }

}

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

//__________________________________________________________________________________________________
void MvfConversions::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 MvfConversions::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 ( c->Charge() != 0 )
          failed.SetBit(eCharge);
                   
        Bool_t failedElectrons=false;
        for (UInt_t i=0; i<c->NDaughters(); ++i) {
          const DaughterData* d = c->DaughterDat(i);
          Bool_t trackQuality = ( i==0 || !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 MvfConversions::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 MvfConversions::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(eEMissedHits);

        if (s->NWrongHits() > fWrongHitsMax)
          failed.SetBit(eEWrongHits);
        
        if (trackQuality && t->NHits()<fNHitsMin)
          failed.SetBit(eENHits);
  
        if (trackQuality && t->Prob()<fTrackProbMin)
          failed.SetBit(eEProb);
          
        if (trackQuality && fExcludePXB1 && t->Hit(Track::PXB1))
          failed.SetBit(eEPxb1);

        return failed;
}

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

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

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

}

BitMask16 MvfConversions::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* MvfConversions::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* MvfConversions::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* MvfConversions::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* MvfConversions::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 MvfConversions::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 MvfConversions::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 MvfConversions::PassPhotonCuts(const Photon *p)
{
  if (p->IsConverted())
    return false;
  
  return true;
}

Double_t MvfConversions::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 MvfConversions::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 MvfConversions::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 *MvfConversions::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 *MvfConversions::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 *MvfConversions::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 MvfConversions::AddRhoRange(Double_t lower, Double_t upper)
{
  fRhoLbs.push_back(lower);
  fRhoUbs.push_back(upper);
}

Bool_t MvfConversions::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;
}
Revision:	1.4
Committed:	Wed Dec 16 18:10:38 2009 UTC (15 years, 4 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.3:	+174 -18 lines
Log Message:	updated MitConversions