Mods/src/ElectronIDMod.cc

// $Id: ElectronIDMod.cc,v 1.24 2009/06/01 17:31:38 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),
  fReverseIsoCut(kFALSE),
  fReverseD0Cut(kFALSE),
  fElIdType(kIdUndef),
  fElIsoType(kIsoUndef),
  fElectrons(0),
  fConversions(0),
  fVertices(0)
{
  // 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;
        
          if     (fReverseIsoCut == kTRUE &&
                  isocut == kFALSE && totalIso < 10)
            isocut = kTRUE;
          else if(fReverseIsoCut == kTRUE)
            isocut = kFALSE;
        }
        break;
      case kNoIso:
        isocut = kTRUE;
        break;
      case kCustomIso:
      default:
        break;
    }

    if (isocut == kFALSE)
      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) {
      Bool_t d0cut = kFALSE;
      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) d0cut = kTRUE;

      if     (fReverseD0Cut == kTRUE &&
              d0cut == kFALSE && d0_real < 0.05)
        d0cut = kTRUE;
      else if(fReverseD0Cut == kTRUE)
        d0cut = kFALSE;

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