Mods/src/MvfConversions.cc

// $Id: FullExample.cc,v 1.2 2008/06/09 11:47:03 paus 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) :
  TAModule    (name,title),
  fPtMin   (0),
  fEtaMax  (0),
  fProbMin (0),
  fIsoRadius(0),
  fIsoMax(0),
  fExcludePXB1(false),
  fMassMax(0),
  fMissedHitsMax(0),
  fWrongHitsMax(0),
  fElectronPtMin(0),
  fElectronEtaMax(0),
  fComputeEff(kFALSE),
  fTracks     (0),
  fGsfTracks  (0),
  fInOutTracks (0),
  fOutInTracks (0),
  fElectrons  (0),
  fPhotons (0),
  fMCParticles (0),
  fMvfConversions (0),
  fMvfConversionsUnconstrained (0),
  fMCPartName (Names::gkMCPartBrn),
  fTrackName  (Names::gkTrackBrn),
  fConvElectronName   ("Electrons"),
  fPhotonName   (Names::gkPhotonBrn),
  fMvfConvName("MvfConversions"),
  fMvfConvUnconstrainedName("MvfConversionsUnconstrained"),
  fElectronName(Names::gkElectronBrn),
  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),
  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(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);
  
  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 (fComputeEff) {
    LoadBranch(fMCPartName);
    for (UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
      MCParticle* p = fMCParticles->At(i);
      UInt_t nFailed = FailedSimCuts(p);
      if (nFailed<=1) {
        if (nFailed==0 || p->Pt() < fPtMin)
          hSimPt->Fill(p->Pt());
        if ( nFailed==0 || !PassRho(p->DecayVertex().Rho()) ) {
            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();
          }
        }
        if ( nFailed==0 || minPt<fElectronPtMin )
          hSimMinElectronPt->Fill(minPt);
      }
    }
  }
 
//   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");

  LoadBranch(fMvfConvName);
  ObjArray<DecayParticle> goodConversions;
  ObjArray<Electron> goodElectrons;
  UInt_t nConversions = 0;
  for (UInt_t i=0; i<fMvfConversions->GetEntries(); ++i) {
    DecayParticle* c = fMvfConversions->At(i);
    UInt_t nFailed = FailedCuts(c);
    
    //make N-2 plots for decay length
    if (nFailed<=2) {
      UInt_t nAllowedFails=0;
      if (!PassRho(c->Position().Rho()))
        nAllowedFails++;
      if (  c->Lxy() < fLxyMin || c->Lz() < fLzMin )
        nAllowedFails++;
        
      if ( (nFailed-nAllowedFails)==0 ) {
        hLxy->Fill(c->Lxy());
        hLxyOverLxyErr->Fill(c->Lxy()/c->LxyError());
        hLzOverLzErr->Fill(c->Lz()/c->LzError());      
        hLz->Fill(c->Lz());
      }
    
    }
    
    //make N-1 plots
    if (nFailed <= 1) {
      const MCParticle *simPhoton = SimMatch(c);
      UInt_t nSimFailed=0;
      if (simPhoton)
        nSimFailed = FailedSimCuts(simPhoton);
      if (nSimFailed>1)
        simPhoton=0;
        
    
      if (nFailed==0 || c->Pt()<fPtMin) {
        hGammaPt->Fill(c->Pt());
        if (simPhoton)
          if (nSimFailed==0 || simPhoton->Pt()<fPtMin) {
            hSimMatchedGammaPt->Fill(c->Pt());
            hSimMatchedSimGammaPt->Fill(simPhoton->Pt()); 
          }
      }
      if (nFailed==0 || TMath::Abs(c->Eta())>fEtaMax) {
        hGammaEta->Fill(c->Eta());
        if (simPhoton)
          if (nSimFailed==0 || TMath::Abs(simPhoton->Eta())>fEtaMax) {
            hSimMatchedGammaEta->Fill(c->Eta());
            hSimMatchedSimGammaEta->Fill(simPhoton->Eta());
          }
      }
      if ( nFailed==0 || !PassRho(c->Position().Rho()) ) {
        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)
          if ( nSimFailed==0 || !PassRho(simPhoton->DecayVertex().Rho()) ) {
            hSimMatchedConvRadius->Fill(c->Position().Rho());
            hSimMatchedSimConvRadius->Fill(simPhoton->DecayVertex().Rho());
          }
            
      }
      if (nFailed==0 || c->Prob() < fProbMin) {
        hConvProb->Fill(c->Prob());
        hConvChi2->Fill(c->Chi2());
          if (simPhoton)
            if( nSimFailed==0)
              hSimMatchedConvProb->Fill(c->Prob());
      }
      if ( nFailed==0 ) {
        goodConversions.Add(c);
        ++nConversions;
        //fill photon histograms
        hConversionZ->Fill(c->Position().Z());
        hConvDCotTheta->Fill(DCotTheta(c));
        hGammaMass->Fill(c->Mass());
        hDxy->Fill(c->Dxy());
        LoadBranch(fTrackName);
        hNTracks->Fill(NTracks(fTracks));
        Double_t isolation = TrackPtIsolation(c,fIsoRadius);
        hIsolation->Fill(isolation);
        Double_t minElectronPt = c->Daughter(0)->Pt();
        const Particle *minPtElectron = c->Daughter(0);
        for (UInt_t j=0; j<c->NDaughters(); ++j) {
          const Particle* d = c->Daughter(j);
          hElectronPt->Fill(d->Pt());
          hElectronEta->Fill(d->Eta());
          if (d->Pt() < minElectronPt) {
            minElectronPt = d->Pt();
            minPtElectron = d;
          }
          const ChargedParticle* dc = dynamic_cast<const ChargedParticle*>(d);
          if (dc) {
            const Track* t = dc->Trk();
            hTrackChi2->Fill(t->Chi2()/((Double_t)t->Ndof()));
            hTrackProb->Fill(t->Prob());
            hTrackNHits->Fill(t->NHits());
            hTrackNHitsProb->Fill(t->NHits(),t->RChi2());
            hTrackD0->Fill(t->D0());
//             LoadBranch(fElectronName);
//             Electron *e = MatchingElectron(dc,fElectrons);
//             if (e) {
//               //printf("Matching Electron added\n");
//               goodElectrons.Add(e);
//             }
//             else
//               ;//printf("No matching electron found\n");
          }
        }
        hMinElectronPt->Fill(minElectronPt);
        
        if (simPhoton && nSimFailed==0) {

          //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());
          hSimMatchedLxy->Fill(c->Lxy());
          hSimMatchedLz->Fill(c->Lz());
          hSimMatchedDxy->Fill(c->Dxy());
          hSimMatchedNTracks->Fill(NTracks(fTracks));
          hSimMatchedIsolation->Fill(isolation);
          hSimMatchedMinElectronPt->Fill(minElectronPt);
          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();
              hSimMatchedTrackChi2->Fill(t->Chi2()/((Double_t)t->Ndof()));
              hSimMatchedTrackProb->Fill(t->Prob());
              hSimMatchedTrackNHits->Fill(t->NHits());
              hSimMatchedTrackNHitsProb->Fill(t->NHits(),t->RChi2());
              hSimMatchedTrackD0->Fill(t->D0());
              if (d==minPtElectron)
                hSimMatchedMinElectronSimPt->Fill(t->MCPart()->Pt());
            }
          }
        }  
      }
    }

//       else if (!simPhoton) {
//         //hConvDCotThetaPreCut->Fill(DCotTheta(c));
//       // Bool_t toFill = true;
// //        if (TMath::Abs(DCotTheta(c))<0.01)
// //           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 MCParticle *dcmc = dc->Trk()->MCPart();
// //               if (dcmc) {
// //                 ;//printf("Daughter pdg = %i, mother pdg = %i\n",dcmc->PdgId(),dcmc->DistinctMother()->PdgId());
// //                 //if (dcmc->DistinctMother()->PdgId()==22 || TMath::Abs(dcmc->DistinctMother()->PdgId())==3122)
// //                   //toFill=false;
// //               }
// //             }
// //           }
//           //if (toFill)
//             hConvDCotThetaPreCut->Fill(DCotTheta(c));
//             hSimMatchedEPairDeltaPhi->Fill(DPhi(c));
//             hSimMatchedConvChi2->Fill(c->Chi2()/c->Ndof());
//       }
//       else {
//         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();
//             hUnMatchedTrackChi2->Fill(t->Chi2());
//             hUnMatchedTrackProb->Fill(t->Prob());
//             hUnMatchedTrackNHits->Fill(t->NHits());
//             hUnMatchedTrackNHitsProb->Fill(t->NHits(),t->Prob());
//           }
//         }
//       
//       }
//    }
    
  
  }
  hNConversions->Fill(nConversions);
  
//   LoadBranch(fPhotonName);
  for (UInt_t i=0; i<goodConversions.GetEntries(); i++) {
    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.
}

UInt_t MvfConversions::FailedCuts(const DecayParticle* c) {


        //if (p->GetVertex()->Prob() < 0.005)
        //      return false;
        UInt_t nFailed = 0;
        
        if ( c->NDaughters() != 2)
          nFailed++;
          
        if ( c->Pt() < fPtMin )
          nFailed++;
          
        if ( !PassRho(c->Position().Rho()) )
          nFailed++;
          
        if ( TMath::Abs(c->Eta()) > fEtaMax )
          nFailed++;
          
          
        if ( c->Lxy() < fLxyMin || c->Lz() < fLzMin )
          nFailed++;
          
        if ( (c->Lxy()/c->LxyError()) < fLxyOverLxyErrMin )
          nFailed++;
          
        if ( (c->Lz()/c->LzError()) < fLzOverLzErrMin )
          nFailed++;
          
        if ( TMath::Abs(c->Dxy()) > fAbsDxyMax )
          nFailed++;

        if (c->Mass() > fMassMax)
          nFailed++;
          
        if (nFailed>=3)
          return nFailed;
          
        if ( c->Charge() != 0 )
          nFailed++;
          
        if (nFailed>=3)
          return nFailed;
          
        Bool_t failedElectrons=false;
        for (UInt_t i=0; i<c->NDaughters(); ++i) {
          const DaughterData* d = c->DaughterDat(i);
          if (!PassElectronCuts(d))
            failedElectrons=true;
        }
        if (failedElectrons)
          nFailed++;
        
        Double_t isolation = TrackPtIsolation(c,fIsoRadius);
        
        if (isolation > fIsoMax)
          nFailed++;         
        
        if ( c->Prob() < fProbMin )
          nFailed++;

        return nFailed;
//         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().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;
}

Bool_t MvfConversions::PassElectronCuts(const DaughterData* d) {
        //if (p->Pt() < 5.0)
        //  return false;
        
        if (d->Pt() < fElectronPtMin)
          return false;
          
        if (TMath::Abs(d->Eta()) > fElectronEtaMax)
          return false;
        
        const StableData *s = dynamic_cast<const StableData*>(d);

        if (!s)
          return false;

        if (s->NMissedHits() > fMissedHitsMax)
          return false;

        if (s->NWrongHits() > fWrongHitsMax)
          return false;

        const ChargedParticle* c = dynamic_cast<const ChargedParticle*>(s->Original());
        
        if (!c)
          return false;
          
        const Track* t = c->Trk();
        
        if (!t)
          return false;
          
        if (!PassTrackCuts(t))
          return false;
          
        return true;
}

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

 // if ( t->Prob()<fTrackProbMin && t->NHits()<fNHitsMin )
  if (t->NHits()<fNHitsMin)
    return false;
  
  if (t->Prob()<fTrackProbMin)
    return false;  
    
  if (fExcludePXB1 && t->Hit(Track::PXB1))
    return false;

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

}

UInt_t MvfConversions::FailedSimCuts(const MCParticle* p) {
        
        UInt_t nFailed = 0;
        
        if (p->PdgId()!=22)
                nFailed++;
        if (TMath::Abs(p->Eta())>fEtaMax)
                nFailed++;
        //if (p->Pt()<5.0 || p->Pt()>25.0)
        //if (p->Pt()<25.0)
        if (p->Pt()<fPtMin)
                nFailed++;
                
        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) )
          nFailed++;
                
        if (nFailed>=2)
          return nFailed;
          
        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)
          nFailed++;
        if (failedElectronEta)
          nFailed++;
        
        if (nEPlus!=1 || nEMinus!=1)
                nFailed++;
        
        return nFailed;
}

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 = electron1->Trk()->MCPart();
  const MCParticle *simElectron2 = electron2->Trk()->MCPart();
  
  
  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;
  

}

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;

  ThreeVector mom1 =  ((const StableData*)c->DaughterDat(0))->ThreeMom();
  ThreeVector 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;
}

Electron *MvfConversions::MatchingElectron(const ChargedParticle *c, ElectronCol *col) const
{
   Electron *match=0;
   Double_t dRMin=0;
   for (UInt_t i=0; i<col->GetEntries(); ++i) {
     Electron *e = col->At(i);
     Double_t dR = MathUtils::DeltaR(c->Mom(), e->Mom());
     if ( dR<0.1 && (!match || dR<dRMin) ) {
       dRMin = dR;
       match = e;
     }
//     if ( c->TrackerTrk() == e->TrackerTrk() )
//       return e;
   }
  
  return match;
}

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.1
Committed:	Wed Nov 19 17:33:26 2008 UTC (16 years, 5 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
Log Message:	Added MitConversions