ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc
Revision: 1.1
Committed: Wed Jun 1 18:11:52 2011 UTC (13 years, 11 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Log Message:
CiC BaseLine Selection

File Contents

# User Rev Content
1 fabstoec 1.1 // $Id: PhotonCiCMod.cc,v 1.21 2011/05/18 14:01:19 bendavid Exp $
2    
3     #include "MitPhysics/Mods/interface/PhotonCiCMod.h"
4     #include "MitAna/DataTree/interface/PhotonCol.h"
5     #include "MitPhysics/Init/interface/ModNames.h"
6     #include "MitPhysics/Utils/interface/IsolationTools.h"
7     #include "MitPhysics/Utils/interface/PhotonTools.h"
8    
9     using namespace mithep;
10    
11     ClassImp(mithep::PhotonCiCMod)
12    
13     //--------------------------------------------------------------------------------------------------
14     PhotonCiCMod::PhotonCiCMod(const char *name, const char *title) :
15     BaseMod(name,title),
16     fPhotonBranchName (Names::gkPhotonBrn),
17     fGoodPhotonsName (ModNames::gkGoodPhotonsName),
18     fTrackBranchName (Names::gkTrackBrn),
19     fPileUpDenName (Names::gkPileupEnergyDensityBrn),
20     fElectronName ("Electrons"),
21     fPhotonPtMin(15.0),
22     fApplySpikeRemoval(kFALSE),
23     fAbsEtaMax(999.99),
24     fPhotons(0),
25     fTracks(0),
26     fPileUpDen(0),
27     fElectrons(0),
28    
29     fPVName("PrimaryVertexes"),
30     fPV(0),
31     fPVFromBranch(kTRUE)
32    
33     {
34     // Constructor.
35     }
36    
37     //--------------------------------------------------------------------------------------------------
38     void PhotonCiCMod::Process()
39     {
40     // Process entries of the tree.
41    
42     LoadEventObject(fPhotonBranchName, fPhotons);
43    
44     Double_t _tRho = -1.;
45     if (fPhotons->GetEntries()>0) {
46     LoadEventObject(fTrackBranchName, fTracks);
47     LoadEventObject(fPileUpDenName, fPileUpDen);
48     LoadEventObject(fElectronName, fElectrons);
49     LoadEventObject(fPVName, fPV);
50    
51     if(fPileUpDen->GetEntries() > 0)
52     _tRho = (Double_t) fPileUpDen->At(0)->Rho();
53    
54     } else return;
55    
56     PhotonOArr *GoodPhotons = new PhotonOArr;
57     GoodPhotons->SetName(fGoodPhotonsName);
58    
59     std::vector<const Photon*> phCat1;
60     std::vector<const Photon*> phCat2;
61     std::vector<const Photon*> phCat3;
62     std::vector<const Photon*> phCat4;
63    
64    
65     const BaseVertex* theVtx = NULL;
66     if(fPV->GetEntries() > 0)
67     theVtx = fPV->At(0);
68     else return;
69    
70     float cic4_allcuts_temp_sublead[] = {
71     3.77459, 2.18305, 1.76811, 1.30029,
72     11.5519, 3.48306, 3.87394, 1.89038,
73     3.47103, 2.1822, 2.26931, 1.43769,
74     0.0105631, 0.00974116, 0.0282572, 0.0265096,
75     0.0834648, 0.0821447, 0.0648449, 0.0476437,
76     0.94, 0.355935, 0.94, 0.316358,
77     98.9979, 1.94497, 96.2292, 96.2294,
78     1.5, 1.5, 1.5, 1.5 }; // THIS IS ????
79    
80     for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
81     const Photon *ph = fPhotons->At(i);
82    
83     if (ph->Pt() <= fPhotonPtMin) continue;
84     if (ph->AbsEta() >= fAbsEtaMax) continue;
85    
86     Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
87    
88     Bool_t passSpikeRemovalFilter = kTRUE;
89    
90     if (ph->SCluster() && ph->SCluster()->Seed()) {
91     if(ph->SCluster()->Seed()->Energy() > 5.0 &&
92     ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
93     passSpikeRemovalFilter = kFALSE;
94     }
95    
96     if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
97    
98     // compute all relevant observables first
99     double ecalIso = ph->EcalRecHitIsoDr03();
100     double hcalIso = ph->HcalTowerSumEtDr03();
101     double trackIso1 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
102     double trackIso2 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, fPV);
103     double trackIso3 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
104    
105     double combIso1 = ecalIso+hcalIso+trackIso1 - 0.17*_tRho;
106     double combIso2 = ecalIso+hcalIso+trackIso2 - 0.52*_tRho;
107    
108     double tIso1 = (combIso1)*50./ph->Et();
109     double tIso2 = (combIso2)*50./ph->Et();
110     double tIso3 = (trackIso3)*50./ph->Et();
111    
112     double covIEtaIEta =ph->CoviEtaiEta();
113     double HoE = ph->HadOverEm();
114    
115     double R9 = ph->R9();
116    
117     double dRTrack = PhotonTools::ElectronVetoCiC(ph, fElectrons);
118    
119     // stupid hard coded, will update later...
120     if( tIso1 < cic4_allcuts_temp_sublead[0+0*4] &&
121     tIso2 < cic4_allcuts_temp_sublead[0+1*4] &&
122     tIso3 < cic4_allcuts_temp_sublead[0+2*4] &&
123     covIEtaIEta < cic4_allcuts_temp_sublead[0+3*4] &&
124     HoE < cic4_allcuts_temp_sublead[0+4*4] &&
125     R9 > cic4_allcuts_temp_sublead[0+5*4] &&
126     dRTrack*100. > cic4_allcuts_temp_sublead[0+6*4] &&
127     isbarrel ) phCat1.push_back(ph);
128    
129     if( tIso1 < cic4_allcuts_temp_sublead[1+0*4] &&
130     tIso2 < cic4_allcuts_temp_sublead[1+1*4] &&
131     tIso3 < cic4_allcuts_temp_sublead[1+2*4] &&
132     covIEtaIEta < cic4_allcuts_temp_sublead[1+3*4] &&
133     HoE < cic4_allcuts_temp_sublead[1+4*4] &&
134     R9 > cic4_allcuts_temp_sublead[1+5*4] &&
135     dRTrack*100. > cic4_allcuts_temp_sublead[1+6*4] &&
136     isbarrel ) phCat2.push_back(ph);
137    
138     if( tIso1 < cic4_allcuts_temp_sublead[2+0*4] &&
139     tIso2 < cic4_allcuts_temp_sublead[2+1*4] &&
140     tIso3 < cic4_allcuts_temp_sublead[2+2*4] &&
141     covIEtaIEta < cic4_allcuts_temp_sublead[2+3*4] &&
142     HoE < cic4_allcuts_temp_sublead[2+4*4] &&
143     R9 > cic4_allcuts_temp_sublead[2+5*4] &&
144     dRTrack*100. > cic4_allcuts_temp_sublead[2+6*4] ) phCat3.push_back(ph);
145    
146     if( tIso1 < cic4_allcuts_temp_sublead[3+0*4] &&
147     tIso2 < cic4_allcuts_temp_sublead[3+1*4] &&
148     tIso3 < cic4_allcuts_temp_sublead[3+2*4] &&
149     covIEtaIEta < cic4_allcuts_temp_sublead[3+3*4] &&
150     HoE < cic4_allcuts_temp_sublead[3+4*4] &&
151     R9 > cic4_allcuts_temp_sublead[3+5*4] &&
152     dRTrack*100. > cic4_allcuts_temp_sublead[3+6*4] ) phCat4.push_back(ph);
153    
154     }
155    
156     // check which category has two vgalid passes
157     bool foundCat=false;
158     if(phCat1.size() > 1) {
159     // pick the first two guys, assuming they are sorted...
160     GoodPhotons->Add(phCat1[0]);
161     GoodPhotons->Add(phCat1[1]);
162     } else if(phCat2.size() > 1) {
163     GoodPhotons->Add(phCat2[0]);
164     GoodPhotons->Add(phCat2[1]);
165     } else {
166     if(phCat3.size() > 1) {
167     // one must be a Endcap photons... (stupid, but that's how it is in CiC)
168     if(phCat3[0]->SCluster()->AbsEta() > 1.5 || phCat3[1]->SCluster()->AbsEta()) {
169     GoodPhotons->Add(phCat3[0]);
170     GoodPhotons->Add(phCat3[1]);
171     foundCat=true;
172     }
173     if (!foundCat) {
174     if(phCat4.size() > 1) {
175     if(phCat4[0]->SCluster()->AbsEta() > 1.5 || phCat4[1]->SCluster()->AbsEta()) {
176     GoodPhotons->Add(phCat4[0]);
177     GoodPhotons->Add(phCat4[1]);
178     }
179     }
180     }
181     }
182     }
183    
184     // sort according to pt
185     GoodPhotons->Sort();
186    
187     // add to event for other modules to use
188     AddObjThisEvt(GoodPhotons);
189     }
190    
191     //--------------------------------------------------------------------------------------------------
192     void PhotonCiCMod::SlaveBegin()
193     {
194     // Run startup code on the computer (slave) doing the actual analysis. Here,
195     // we just request the photon collection branch.
196    
197     ReqEventObject(fPhotonBranchName, fPhotons, kTRUE);
198     ReqEventObject(fTrackBranchName, fTracks, kTRUE);
199     ReqEventObject(fElectronName, fElectrons, kTRUE);
200     ReqEventObject(fPileUpDenName, fPileUpDen, kTRUE);
201     ReqEventObject(fPVName, fPV, fPVFromBranch);
202    
203     }