ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc
Revision: 1.2
Committed: Thu Jun 2 14:18:07 2011 UTC (13 years, 11 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.1: +20 -8 lines
Log Message:
add fiducial cut and a few fixes

File Contents

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