Mods/src/ElectronIDMod.cc

// $Id: ElectronIDMod.cc,v 1.22 2009/04/30 08:10:11 ceballos Exp $

#include "MitPhysics/Mods/interface/ElectronIDMod.h"
#include "MitPhysics/Init/interface/ModNames.h"

using namespace mithep;

ClassImp(mithep::ElectronIDMod)

//--------------------------------------------------------------------------------------------------
ElectronIDMod::ElectronIDMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fElectronBranchName(Names::gkElectronBrn),
  fConversionBranchName(Names::gkMvfConversionBrn),
  fGoodElectronsName(ModNames::gkGoodElectronsName),  
  fVertexName(string("PrimaryVertexesBeamSpot").c_str()),
  fElectronIDType("Tight"),
  fElectronIsoType("TrackJuraSliding"),
  fElectronPtMin(10),
  fIDLikelihoodCut(0.9),
  fTrackIsolationCut(5.0),
  fCaloIsolationCut(5.0),
  fEcalJuraIsoCut(5.0),
  fHcalIsolationCut(5.0),
  fApplyConvFilter(kTRUE),
  fApplyD0Cut(kTRUE),
  fD0Cut(0.025),
  fElIdType(kIdUndef),
  fElIsoType(kIsoUndef),
  fElectrons(0),
  fConversions(0),
  fVertices(0),
  fReverseIsoCut(kFALSE)
{
  // Constructor.
}

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

  LoadEventObject(fElectronBranchName, fElectrons);

  ElectronOArr *GoodElectrons = new ElectronOArr;
  GoodElectrons->SetName(fGoodElectronsName);

  for (UInt_t i=0; i<fElectrons->GetEntries(); ++i) {
    const Electron *e = fElectrons->At(i);        

    if (e->Pt() <= fElectronPtMin) 
      continue;
    
    Bool_t idcut = kFALSE;
    switch (fElIdType) {
      case kTight:
        idcut = e->PassTightID();
        break;
      case kLoose:
        idcut = e->PassLooseID();
       break;
      case kLikelihood:
        idcut = (e->IDLikelihood() > fIDLikelihoodCut);
        break;
      case kNoId:
        idcut = kTRUE;
        break;
      case kCustomId:
      default:
        break;
    }

    if (!idcut) 
      continue;

    Bool_t isocut = kFALSE;
    switch (fElIsoType) {
      case kTrackCalo:
        isocut = (e->TrackIsolation() < fTrackIsolationCut) &&
                 (e->CaloIsolation() < fCaloIsolationCut);
        break;
      case kTrackJura:
        isocut = (e->TrackIsolation() < fTrackIsolationCut) &&
                 (e->EcalJurassicIsolation() < fEcalJuraIsoCut) &&
                 (e->HcalIsolation() < fHcalIsolationCut);
        break;
      case kTrackJuraSliding:
        {
          Double_t totalIso = e->TrackIsolation() + e->EcalJurassicIsolation() - 1.5;
          if ((totalIso < (e->Pt()-10.0)*5.0/15.0 && e->Pt() <= 25) ||
              (totalIso < 5.0 && e->Pt() > 25) ||
               totalIso <= 0)
            isocut = kTRUE;
        }
        break;
      case kNoIso:
        isocut = kTRUE;
        break;
      case kCustomIso:
      default:
        break;
    }

    if ((isocut == kFALSE && fReverseIsoCut == kFALSE) ||
        (isocut == kTRUE  && fReverseIsoCut == kTRUE))
      continue;

    // apply conversion filter
    Bool_t isGoodConversion = kFALSE;
    if (fApplyConvFilter) {
      LoadBranch(fConversionBranchName);
      for (UInt_t ifc=0; ifc<fConversions->GetEntries(); ifc++) {
        
        Bool_t ConversionMatchFound = kFALSE;
        for (UInt_t d=0; d<fConversions->At(ifc)->NDaughters(); d++) {
          const Track *trk = dynamic_cast<const ChargedParticle*>
                       (fConversions->At(ifc)->Daughter(d))->Trk();
          if (e->GsfTrk() == trk) {
            ConversionMatchFound = kTRUE;
            break;
          }
        }

        // if match between the e-track and one of the conversion legs
        if (ConversionMatchFound == kTRUE){
          isGoodConversion =  (fConversions->At(ifc)->Prob() > 0.0005) &&
                              (fConversions->At(ifc)->Lxy() > 0) &&
                              (fConversions->At(ifc)->Lz() > 0);

          if (isGoodConversion == kTRUE) {
            for (UInt_t d=0; d<fConversions->At(ifc)->NDaughters(); d++) {
              const Track *trk = dynamic_cast<const ChargedParticle*>
                           (fConversions->At(ifc)->Daughter(d))->Trk();
                      
              if (trk) {
                // These requirements are not used for the GSF track (d == 1)
                if (!(trk->NHits() > 8 && trk->Prob() > 0.005) && d == 0)
                  isGoodConversion = kFALSE;
              
                const StableData *sd = dynamic_cast<const StableData*>
                                  (fConversions->At(ifc)->DaughterDat(d));
                if (sd->NWrongHits() != 0)
                  isGoodConversion = kFALSE;
              
              } else {
                isGoodConversion = kFALSE;
              }
            }
          }
        }

        if (isGoodConversion == kTRUE) break;

      } // loop over all conversions 
      
    }
    if (isGoodConversion == kTRUE) continue;

    if (fApplyD0Cut) {
      LoadBranch(fVertexName);
      // d0 cut
      double d0_real = 99999;
      for(uint i0 = 0; i0 < fVertices->GetEntries(); i0++) {
        double pD0 = e->GsfTrk()->D0Corrected(*fVertices->At(i0));
        if(TMath::Abs(pD0) < TMath::Abs(d0_real)) d0_real = TMath::Abs(pD0);
      }
      if(d0_real >= fD0Cut) continue;
    }

    // add good electron
    GoodElectrons->Add(e);
  }

  // sort according to pt
  GoodElectrons->Sort();

  // add to event for other modules to use
  AddObjThisEvt(GoodElectrons);  
}

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

  ReqEventObject(fElectronBranchName, fElectrons, kTRUE);

  if (fApplyConvFilter)
    ReqEventObject(fConversionBranchName, fConversions, kTRUE);

  if (fApplyD0Cut)
    ReqEventObject(fVertexName, fVertices, kTRUE);

  if (fElectronIDType.CompareTo("Tight") == 0) 
    fElIdType = kTight;
  else if (fElectronIDType.CompareTo("Loose") == 0) 
    fElIdType = kLoose;
  else if (fElectronIDType.CompareTo("Likelihood") == 0) 
    fElIdType = kLikelihood;
  else if (fElectronIDType.CompareTo("NoId") == 0) 
    fElIdType = kNoId;
  else if (fElectronIDType.CompareTo("Custom") == 0) {
    fElIdType = kCustomId;
    SendError(kWarning, "SlaveBegin",
              "Custom electron identification is not yet implemented.");
  } else {
    SendError(kAbortAnalysis, "SlaveBegin",
              "The specified electron identification %s is not defined.",
              fElectronIDType.Data());
    return;
  }

  if (fElectronIsoType.CompareTo("TrackCalo") == 0 )
    fElIsoType = kTrackCalo;
  else if (fElectronIsoType.CompareTo("TrackJura") == 0) 
    fElIsoType = kTrackJura;
  else if(fElectronIsoType.CompareTo("TrackJuraSliding") == 0)
    fElIsoType = kTrackJuraSliding;
  else if (fElectronIsoType.CompareTo("NoIso") == 0 )
    fElIsoType = kNoIso;
  else if (fElectronIsoType.CompareTo("Custom") == 0 ) {
    fElIsoType = kCustomIso;
    SendError(kWarning, "SlaveBegin",
              "Custom electron isolation is not yet implemented.");
  } else {
    SendError(kAbortAnalysis, "SlaveBegin",
              "The specified electron isolation %s is not defined.",
              fElectronIsoType.Data());
    return;
  }
}
Revision:	1.23
Committed:	Tue May 19 11:42:20 2009 UTC (15 years, 11 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
Changes since 1.22:	+5 -5 lines
Log Message:	Use the new ReqEventObject model.
#	User	Rev	Content
1	loizides	1.23	// $Id: ElectronIDMod.cc,v 1.22 2009/04/30 08:10:11 ceballos Exp $
2	loizides	1.1
3			#include "MitPhysics/Mods/interface/ElectronIDMod.h"
4	loizides	1.5	#include "MitPhysics/Init/interface/ModNames.h"
5	loizides	1.1
6			using namespace mithep;
7
8			ClassImp(mithep::ElectronIDMod)
9
10			//--------------------------------------------------------------------------------------------------
11	loizides	1.5	ElectronIDMod::ElectronIDMod(const char name, const char title) :
12	loizides	1.1	BaseMod(name,title),
13	loizides	1.5	fElectronBranchName(Names::gkElectronBrn),
14	ceballos	1.12	fConversionBranchName(Names::gkMvfConversionBrn),
15	loizides	1.5	fGoodElectronsName(ModNames::gkGoodElectronsName),
16	ceballos	1.12	fVertexName(string("PrimaryVertexesBeamSpot").c_str()),
17	loizides	1.1	fElectronIDType("Tight"),
18	ceballos	1.8	fElectronIsoType("TrackJuraSliding"),
19	loizides	1.1	fElectronPtMin(10),
20			fIDLikelihoodCut(0.9),
21			fTrackIsolationCut(5.0),
22			fCaloIsolationCut(5.0),
23	loizides	1.5	fEcalJuraIsoCut(5.0),
24			fHcalIsolationCut(5.0),
25	loizides	1.14	fApplyConvFilter(kTRUE),
26	ceballos	1.18	fApplyD0Cut(kTRUE),
27	loizides	1.14	fD0Cut(0.025),
28	loizides	1.5	fElIdType(kIdUndef),
29	ceballos	1.12	fElIsoType(kIsoUndef),
30	loizides	1.14	fElectrons(0),
31			fConversions(0),
32	ceballos	1.21	fVertices(0),
33	ceballos	1.22	fReverseIsoCut(kFALSE)
34	loizides	1.1	{
35			// Constructor.
36			}
37
38			//--------------------------------------------------------------------------------------------------
39			void ElectronIDMod::Process()
40			{
41			// Process entries of the tree.
42
43	loizides	1.23	LoadEventObject(fElectronBranchName, fElectrons);
44	loizides	1.1
45	loizides	1.6	ElectronOArr *GoodElectrons = new ElectronOArr;
46			GoodElectrons->SetName(fGoodElectronsName);
47	loizides	1.1
48	ceballos	1.18	for (UInt_t i=0; i<fElectrons->GetEntries(); ++i) {
49	loizides	1.5	const Electron *e = fElectrons->At(i);
50	loizides	1.1
51	loizides	1.5	if (e->Pt() <= fElectronPtMin)
52			continue;
53	loizides	1.1
54	loizides	1.5	Bool_t idcut = kFALSE;
55			switch (fElIdType) {
56			case kTight:
57			idcut = e->PassTightID();
58			break;
59			case kLoose:
60			idcut = e->PassLooseID();
61			break;
62			case kLikelihood:
63			idcut = (e->IDLikelihood() > fIDLikelihoodCut);
64			break;
65	loizides	1.11	case kNoId:
66			idcut = kTRUE;
67			break;
68	loizides	1.5	case kCustomId:
69			default:
70			break;
71	loizides	1.1	}
72
73	loizides	1.5	if (!idcut)
74			continue;
75
76			Bool_t isocut = kFALSE;
77			switch (fElIsoType) {
78			case kTrackCalo:
79			isocut = (e->TrackIsolation() < fTrackIsolationCut) &&
80			(e->CaloIsolation() < fCaloIsolationCut);
81			break;
82			case kTrackJura:
83			isocut = (e->TrackIsolation() < fTrackIsolationCut) &&
84			(e->EcalJurassicIsolation() < fEcalJuraIsoCut) &&
85			(e->HcalIsolation() < fHcalIsolationCut);
86			break;
87			case kTrackJuraSliding:
88	ceballos	1.18	{
89	loizides	1.5	Double_t totalIso = e->TrackIsolation() + e->EcalJurassicIsolation() - 1.5;
90	ceballos	1.12	if ((totalIso < (e->Pt()-10.0)*5.0/15.0 && e->Pt() <= 25) \|\|
91			(totalIso < 5.0 && e->Pt() > 25) \|\|
92	ceballos	1.7	totalIso <= 0)
93	loizides	1.5	isocut = kTRUE;
94			}
95			break;
96			case kNoIso:
97			isocut = kTRUE;
98			break;
99			case kCustomIso:
100			default:
101			break;
102	loizides	1.1	}
103
104	ceballos	1.21	if ((isocut == kFALSE && fReverseIsoCut == kFALSE) \|\|
105			(isocut == kTRUE && fReverseIsoCut == kTRUE))
106	loizides	1.5	continue;
107
108	loizides	1.14	// apply conversion filter
109	ceballos	1.12	Bool_t isGoodConversion = kFALSE;
110	loizides	1.14	if (fApplyConvFilter) {
111	ceballos	1.12	LoadBranch(fConversionBranchName);
112			for (UInt_t ifc=0; ifc<fConversions->GetEntries(); ifc++) {
113
114			Bool_t ConversionMatchFound = kFALSE;
115			for (UInt_t d=0; d<fConversions->At(ifc)->NDaughters(); d++) {
116			const Track trk = dynamic_cast<const ChargedParticle>
117			(fConversions->At(ifc)->Daughter(d))->Trk();
118	ceballos	1.13	if (e->GsfTrk() == trk) {
119	ceballos	1.12	ConversionMatchFound = kTRUE;
120			break;
121			}
122			}
123
124			// if match between the e-track and one of the conversion legs
125			if (ConversionMatchFound == kTRUE){
126	ceballos	1.19	isGoodConversion = (fConversions->At(ifc)->Prob() > 0.0005) &&
127	ceballos	1.12	(fConversions->At(ifc)->Lxy() > 0) &&
128			(fConversions->At(ifc)->Lz() > 0);
129
130			if (isGoodConversion == kTRUE) {
131			for (UInt_t d=0; d<fConversions->At(ifc)->NDaughters(); d++) {
132			const Track trk = dynamic_cast<const ChargedParticle>
133			(fConversions->At(ifc)->Daughter(d))->Trk();
134
135			if (trk) {
136	ceballos	1.20	// These requirements are not used for the GSF track (d == 1)
137			if (!(trk->NHits() > 8 && trk->Prob() > 0.005) && d == 0)
138	ceballos	1.12	isGoodConversion = kFALSE;
139
140			const StableData sd = dynamic_cast<const StableData>
141			(fConversions->At(ifc)->DaughterDat(d));
142			if (sd->NWrongHits() != 0)
143			isGoodConversion = kFALSE;
144
145			} else {
146			isGoodConversion = kFALSE;
147			}
148			}
149			}
150			}
151
152			if (isGoodConversion == kTRUE) break;
153
154			} // loop over all conversions
155
156			}
157			if (isGoodConversion == kTRUE) continue;
158
159	ceballos	1.15	if (fApplyD0Cut) {
160			LoadBranch(fVertexName);
161			// d0 cut
162			double d0_real = 99999;
163			for(uint i0 = 0; i0 < fVertices->GetEntries(); i0++) {
164			double pD0 = e->GsfTrk()->D0Corrected(*fVertices->At(i0));
165			if(TMath::Abs(pD0) < TMath::Abs(d0_real)) d0_real = TMath::Abs(pD0);
166			}
167			if(d0_real >= fD0Cut) continue;
168	ceballos	1.12	}
169
170	loizides	1.5	// add good electron
171	ceballos	1.12	GoodElectrons->Add(e);
172	loizides	1.5	}
173	loizides	1.1
174	loizides	1.9	// sort according to pt
175			GoodElectrons->Sort();
176
177	loizides	1.5	// add to event for other modules to use
178	loizides	1.6	AddObjThisEvt(GoodElectrons);
179	loizides	1.1	}
180
181			//--------------------------------------------------------------------------------------------------
182			void ElectronIDMod::SlaveBegin()
183			{
184			// Run startup code on the computer (slave) doing the actual analysis. Here,
185	loizides	1.5	// we just request the electron collection branch.
186	loizides	1.1
187	loizides	1.23	ReqEventObject(fElectronBranchName, fElectrons, kTRUE);
188	loizides	1.17
189	ceballos	1.16	if (fApplyConvFilter)
190	loizides	1.23	ReqEventObject(fConversionBranchName, fConversions, kTRUE);
191	loizides	1.17
192	ceballos	1.15	if (fApplyD0Cut)
193	loizides	1.23	ReqEventObject(fVertexName, fVertices, kTRUE);
194	loizides	1.1
195	loizides	1.5	if (fElectronIDType.CompareTo("Tight") == 0)
196			fElIdType = kTight;
197			else if (fElectronIDType.CompareTo("Loose") == 0)
198			fElIdType = kLoose;
199			else if (fElectronIDType.CompareTo("Likelihood") == 0)
200			fElIdType = kLikelihood;
201	loizides	1.10	else if (fElectronIDType.CompareTo("NoId") == 0)
202			fElIdType = kNoId;
203	loizides	1.5	else if (fElectronIDType.CompareTo("Custom") == 0) {
204			fElIdType = kCustomId;
205			SendError(kWarning, "SlaveBegin",
206			"Custom electron identification is not yet implemented.");
207			} else {
208			SendError(kAbortAnalysis, "SlaveBegin",
209			"The specified electron identification %s is not defined.",
210			fElectronIDType.Data());
211			return;
212			}
213
214			if (fElectronIsoType.CompareTo("TrackCalo") == 0 )
215			fElIsoType = kTrackCalo;
216			else if (fElectronIsoType.CompareTo("TrackJura") == 0)
217			fElIsoType = kTrackJura;
218			else if(fElectronIsoType.CompareTo("TrackJuraSliding") == 0)
219			fElIsoType = kTrackJuraSliding;
220			else if (fElectronIsoType.CompareTo("NoIso") == 0 )
221			fElIsoType = kNoIso;
222			else if (fElectronIsoType.CompareTo("Custom") == 0 ) {
223			fElIsoType = kCustomIso;
224			SendError(kWarning, "SlaveBegin",
225			"Custom electron isolation is not yet implemented.");
226			} else {
227			SendError(kAbortAnalysis, "SlaveBegin",
228			"The specified electron isolation %s is not defined.",
229			fElectronIsoType.Data());
230			return;
231			}
232	loizides	1.1	}