ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/rootEWKanalyzer/src/baseClass.C
Revision: 1.6
Committed: Tue Dec 21 15:25:39 2010 UTC (14 years, 4 months ago) by jueugste
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.5: +19 -13 lines
Log Message: 
Added generated PU info

File Contents

# User Rev Content
1 jueugste 1.1 #define baseClass_cxx
2     #include "baseClass.h"
3    
4 jueugste 1.2 #include "TH1I.h"
5     #include "TFile.h"
6    
7 jueugste 1.1 baseClass::baseClass(string * inputList, string * cutFile, string * treeName, string * outputFileName, string * cutEfficFile)
8     {
9     std::cout << "baseClass::baseClass(): begins " << std::endl;
10     inputList_ = inputList;
11     cutFile_ = cutFile;
12     treeName_= treeName;
13     outputFileName_ = outputFileName;
14     cutEfficFile_ = cutEfficFile;
15     init();
16     std::cout << "baseClass::baseClass(): ends " << std::endl;
17     }
18    
19     baseClass::~baseClass()
20     {
21     std::cout << "baseClass::~baseClass(): begin " << std::endl;
22 jueugste 1.2 // the following lines are from the jet prompt thing...
23     // if( !writeCutHistos() )
24     // {
25     // STDOUT("ERROR: writeCutHistos did not complete successfully.");
26     // }
27     // if( !writeCutEfficFile() )
28     // {
29     // STDOUT("ERROR: writeStatFile did not complete successfully.");
30     // }
31 jueugste 1.1 output_root_->Close();
32     std::cout << "baseClass::~baseClass(): end " << std::endl;
33     }
34    
35    
36     double baseClass::multiply(double a, double b)
37     {
38     double c;
39     c = a * b;
40     return c;
41     }
42    
43 jueugste 1.3 void baseClass::scaleHisto(TH1D* histo, double factor) {
44     int binnumber=histo->GetNbinsX();
45     for(int i=0;i<=binnumber+1;i++) {
46     double bincontent=histo->GetBinContent(i);
47     double newbincontent=bincontent*factor;
48     histo->SetBinContent(i,newbincontent);
49     double binerror=histo->GetBinError(i);
50     double newbinerror=binerror*factor;
51     histo->SetBinError(i,newbinerror);
52     }
53     }
54    
55    
56    
57     string baseClass::getFileName() {
58     if(tree_ == NULL){
59     std::cout << "baseClass::init(): ERROR: tree_ = NULL " << std::endl;
60     // return 'NIL';
61     }
62    
63     tree_->GetEntry();
64     TFile *f = tree_->GetCurrentFile();
65     // cout<<"-----------------------------"<<endl;
66     // cout<<f->GetName()<<endl;
67     string filename = f->GetName();
68 jueugste 1.4 cout<<filename<<endl;
69 jueugste 1.3 return filename;
70     }
71    
72    
73 jueugste 1.2 void baseClass::getHLTtable() {
74     bool printTriggerTable=true;
75     if(tree_ == NULL){
76     std::cout << "baseClass::init(): ERROR: tree_ = NULL " << std::endl;
77     return;
78     }
79    
80     tree_->GetEntry();
81     TFile *f = tree_->GetCurrentFile();
82    
83 jueugste 1.3 cout<<f->GetName()<<endl;
84 jueugste 1.2 TH1I *hlt_stats = (TH1I*)f->Get("analyze/HLTTriggerStats");
85    
86     // clear the map
87     if(!HLTmap.empty()) {
88     HLTmap.erase(HLTmap.begin(),HLTmap.end());
89     }
90    
91     for( int i=0; i < hlt_stats->GetNbinsX(); i++ ){
92     if(printTriggerTable==true) cout<<i<<": "<<hlt_stats->GetXaxis()->GetBinLabel(i+1)<<endl;
93     HLTmap[hlt_stats->GetXaxis()->GetBinLabel(i+1)] = i;
94     }
95     }
96    
97 jueugste 1.5 void baseClass::GetHLTNames(Int_t& run){
98    
99     TFile *f = tree_->GetCurrentFile();
100     TTree* runTree = (TTree*)f->Get("analyze/RunInfo");
101     cout<<"here"<<endl;
102     std::vector<std::string>* HLTNames;
103     if ( !runTree ) {
104     std::cerr << "!!! UserAnalysisBase::GetHLTNames "
105     << "Couldn't get analyze/RunInfo tree" << std::endl;
106     return;
107     }
108     //TH1I *hlt_stats = (TH1I*)f->Get("analyze/HLTTriggerStats");
109    
110     /*if ( fVerbose>0 ) */std::cout << "Retrieving HLTNames for run " << run << std::endl;
111     runTree->SetBranchAddress("HLTNames",&HLTNames);
112     runTree->GetEntryWithIndex(run);
113    
114     cout<<HLTNames->size()<<endl;
115    
116     for( int i=0; i < HLTNames->size(); i++ ) {
117     HLTmap[(*HLTNames)[i]] = i;
118     cout<<i<<": "<<(*HLTNames)[i]<<endl;
119     }
120     }
121    
122    
123 jueugste 1.2 int baseClass::getHLTtriggerBit(string triggerName) {
124     if(HLTmap.empty()) {
125 jueugste 1.5 bool verbose = true;
126     if(verbose) cerr<<"HLTtrigger: no HLTmap available!"<<endl;
127 jueugste 1.2 return -1;
128     }
129 jueugste 1.5
130 jueugste 1.2 map<string, int>::iterator iter = HLTmap.find(triggerName);
131     //cout<<"HLT bit: "<<iter->second<<endl; // iter->second is the HLT trigger bit
132     if(iter==HLTmap.end()) {
133 jueugste 1.5 //cerr<<"HLTtrigger: trigger name "<<triggerName<< " not found!"<<endl;
134 jueugste 1.2 return -1;
135     }
136     else {
137     return iter->second;
138     }
139     }
140    
141 jueugste 1.1 bool baseClass::triggerBitSelection(int hlt[200], int l1phys[128], int l1tech[64])
142     {
143     bool passed=false;
144     if((l1tech[40]==1 || l1tech[41]==1) && // pass BSC
145     (l1tech[36]==0 && l1tech[37]==0 && l1tech[38]==0 && l1tech[39]==0)// && // Beam Halo Veto
146     // l1tech[0]==1 && // pass BPTX
147     // hlt[116]==1 // pass Physics Declared Bit
148     )
149     {
150     passed=true;
151     }
152     return passed;
153    
154     }
155    
156     bool baseClass::goodPrimaryVertexSelection(int ndof, double z)
157     {
158     bool passed=false;
159     if(ndof>=5 && fabs(z)<15) passed=true;
160     return passed;
161     }
162    
163     // bool baseClass::noBeamScrapingSelection(int nTracks, int goodTrack)
164     // {
165     // passed=false;
166     // if(nTracks<10) passed=true;
167     // else if(ntracks>=10) {
168     // for(int i=0;i<nTracks;i++) {
169     // if()
170     // }
171     // }
172     // return passed;
173     // }
174    
175     bool myEvent::WorkingPointID(myEvent electron, vector<double> cuts)
176     {
177     bool passed=false;
178     double etaMax_B=1.44;
179     double etaMax_E=2.6;
180     // get the variables
181     double See=electron.sigmaIetaIeta[0];
182     double Dphi=fabs(electron.deltaPhi[0]);
183     double Deta=fabs(electron.deltaEta[0]);
184     double HoE=electron.HoverE[0];
185     double TrkIso=electron.trkIso[0];
186     double EcalIso=electron.ecalIso[0];
187     double HcalIso=electron.hcalIso[0];
188     double missingHits=electron.numberOfMissingInnerHits[0];
189     // get the cuts
190     double See_B=cuts[0];
191     double Dphi_B=cuts[1];
192     double Deta_B=cuts[2];
193     double HoE_B=cuts[3];
194     double TrkIso_B=cuts[4];
195     double EcalIso_B=cuts[5];
196     double HcalIso_B=cuts[6];
197     double missingHits_B=cuts[7];
198     double See_E=cuts[8];
199     double Dphi_E=cuts[9];
200     double Deta_E=cuts[10];
201     double HoE_E=cuts[11];
202     double TrkIso_E=cuts[12];
203     double EcalIso_E=cuts[13];
204     double HcalIso_E=cuts[14];
205     double missingHits_E=cuts[15];
206     if(fabs(electron.fourvector[0].Eta())<etaMax_B) {
207     if(See<See_B &&
208     Dphi<Dphi_B &&
209     Deta<Deta_B &&
210     HoE<HoE_B &&
211     TrkIso<TrkIso_B &&
212     EcalIso<EcalIso_B &&
213     HcalIso<HcalIso_B &&
214     missingHits<missingHits_B) passed=true;
215     }
216     else if(fabs(electron.fourvector[0].Eta())>etaMax_B && fabs(electron.fourvector[0].Eta())<etaMax_E) {
217     if(See<See_E &&
218     Dphi<Dphi_E &&
219     Deta<Deta_E &&
220     HoE<HoE_E &&
221     TrkIso<TrkIso_E &&
222     EcalIso<EcalIso_E &&
223     HcalIso<HcalIso_E &&
224     missingHits<missingHits_E) passed=true;
225     }
226     return passed;
227     }
228    
229     bool myEvent::WorkingPointNminus1(myEvent electron, vector<double> cuts, char *cut)
230     {
231     bool passed=false;
232     int combination[8]={0,0,0,0,0,0,0,0};
233     double etaMax_B=1.44;
234     double etaMax_E=2.6;
235     // get the variables
236     double See=electron.sigmaIetaIeta[0];
237     double Dphi=fabs(electron.deltaPhi[0]);
238     double Deta=fabs(electron.deltaEta[0]);
239     double HoE=electron.HoverE[0];
240     double TrkIso=electron.trkIso[0];
241     double EcalIso=electron.ecalIso[0];
242     double HcalIso=electron.hcalIso[0];
243     double missingHits=electron.numberOfMissingInnerHits[0];
244 jueugste 1.2
245    
246 jueugste 1.1 // get the cuts
247     double See_B=cuts[0];
248     double Dphi_B=cuts[1];
249     double Deta_B=cuts[2];
250     double HoE_B=cuts[3];
251     double TrkIso_B=cuts[4];
252     double EcalIso_B=cuts[5];
253     double HcalIso_B=cuts[6];
254 jueugste 1.2 double CombIso_B=cuts[7];
255     double missingHits_B=cuts[16];
256 jueugste 1.1 double See_E=cuts[8];
257     double Dphi_E=cuts[9];
258     double Deta_E=cuts[10];
259     double HoE_E=cuts[11];
260     double TrkIso_E=cuts[12];
261     double EcalIso_E=cuts[13];
262     double HcalIso_E=cuts[14];
263 jueugste 1.2 double CombIso_E=cuts[15];
264     double missingHits_E=cuts[16];
265     double dist=cuts[17];
266    
267    
268    
269    
270    
271 jueugste 1.1 // get the cut which should NOT be made
272     if(!strcmp(cut,"sigmaIetaIeta")) combination[0]=1;
273     else if(!strcmp(cut,"deltaPhi")) combination[1]=1;
274     else if(!strcmp(cut,"deltaEta")) combination[2]=1;
275     else if(!strcmp(cut,"HoverE")) combination[3]=1;
276     else if(!strcmp(cut,"trkIso")) combination[4]=1;
277     else if(!strcmp(cut,"ecalIso")) combination[5]=1;
278     else if(!strcmp(cut,"hcalIso")) combination[6]=1;
279     else if(!strcmp(cut,"conversions")) combination[7]=1;
280     else cout<<"ERROR: not able to give N-1 flag for the working point!"<<endl;
281     // divide in barrel and endcaps
282     if(fabs(electron.fourvector[0].Eta())<etaMax_B) {
283     // sigma ieta ieta
284     if(combination[0]==1 &&
285     Dphi<Dphi_B &&
286     Deta<Deta_B &&
287     HoE<HoE_B &&
288     TrkIso<TrkIso_B &&
289     EcalIso<EcalIso_B &&
290     HcalIso<HcalIso_B &&
291 jueugste 1.2 missingHits<=missingHits_B) passed=true;
292 jueugste 1.1 // delta phi
293     else if(combination[1]==1 &&
294     See<See_B &&
295     Deta<Deta_B &&
296     HoE<HoE_B &&
297     TrkIso<TrkIso_B &&
298     EcalIso<EcalIso_B &&
299     HcalIso<HcalIso_B &&
300 jueugste 1.2 missingHits<=missingHits_B) passed=true;
301 jueugste 1.1 // delta eta
302     else if(combination[2]==1 &&
303     See<See_B &&
304     Dphi<Dphi_B &&
305     HoE<HoE_B &&
306     TrkIso<TrkIso_B &&
307     EcalIso<EcalIso_B &&
308     HcalIso<HcalIso_B &&
309 jueugste 1.2 missingHits<=missingHits_B) passed=true;
310 jueugste 1.1 // H over E
311     else if(combination[3]==1 &&
312     See<See_B &&
313     Dphi<Dphi_B &&
314     Deta<Deta_B &&
315     TrkIso<TrkIso_B &&
316     EcalIso<EcalIso_B &&
317     HcalIso<HcalIso_B &&
318 jueugste 1.2 missingHits<=missingHits_B) passed=true;
319 jueugste 1.1 // Track isolation
320     else if(combination[4]==1 &&
321     See<See_B &&
322     Dphi<Dphi_B &&
323     Deta<Deta_B &&
324     HoE<HoE_B &&
325     EcalIso<EcalIso_B &&
326     HcalIso<HcalIso_B &&
327 jueugste 1.2 missingHits<=missingHits_B) passed=true;
328 jueugste 1.1 // ECAL isolation
329     else if(combination[5]==1 &&
330     See<See_B &&
331     Dphi<Dphi_B &&
332     Deta<Deta_B &&
333     HoE<HoE_B &&
334     TrkIso<TrkIso_B &&
335     HcalIso<HcalIso_B &&
336 jueugste 1.2 missingHits<=missingHits_B) passed=true;
337 jueugste 1.1 // HCAL isolation
338     else if(combination[6]==1 &&
339     See<See_B &&
340     Dphi<Dphi_B &&
341     Deta<Deta_B &&
342     HoE<HoE_B &&
343     TrkIso<TrkIso_B &&
344     EcalIso<EcalIso_B &&
345 jueugste 1.2 missingHits<=missingHits_B) passed=true;
346 jueugste 1.1 // Conversions
347     else if(combination[7]==1 &&
348     See<See_B &&
349     Dphi<Dphi_B &&
350     Deta<Deta_B &&
351     HoE<HoE_B &&
352     TrkIso<TrkIso_B &&
353     EcalIso<EcalIso_B &&
354     HcalIso<HcalIso_B) passed=true;
355    
356    
357     } // barrel
358     else if(fabs(electron.fourvector[0].Eta())>etaMax_B && fabs(electron.fourvector[0].Eta())<etaMax_E) {
359     // sigma ieta ieta
360     if(combination[0]==1 &&
361     Dphi<Dphi_E &&
362     Deta<Deta_E &&
363     HoE<HoE_E &&
364     TrkIso<TrkIso_E &&
365     EcalIso<EcalIso_E &&
366     HcalIso<HcalIso_E &&
367 jueugste 1.2 missingHits<=missingHits_E) passed=true;
368 jueugste 1.1 // delta phi
369     else if(combination[1]==1 &&
370     See<See_E &&
371     Deta<Deta_E &&
372     HoE<HoE_E &&
373     TrkIso<TrkIso_E &&
374     EcalIso<EcalIso_E &&
375     HcalIso<HcalIso_E &&
376 jueugste 1.2 missingHits<=missingHits_E) passed=true;
377 jueugste 1.1 // delta eta
378     else if(combination[2]==1 &&
379     See<See_E &&
380     Dphi<Dphi_E &&
381     HoE<HoE_E &&
382     TrkIso<TrkIso_E &&
383     EcalIso<EcalIso_E &&
384     HcalIso<HcalIso_E &&
385 jueugste 1.2 missingHits<=missingHits_E) passed=true;
386 jueugste 1.1 // H over E
387     else if(combination[3]==1 &&
388     See<See_E &&
389     Dphi<Dphi_E &&
390     Deta<Deta_E &&
391     TrkIso<TrkIso_E &&
392     EcalIso<EcalIso_E &&
393     HcalIso<HcalIso_E &&
394 jueugste 1.2 missingHits<=missingHits_E) passed=true;
395 jueugste 1.1 // Track isolation
396     else if(combination[4]==1 &&
397     See<See_E &&
398     Dphi<Dphi_E &&
399     Deta<Deta_E &&
400     HoE<HoE_E &&
401     EcalIso<EcalIso_E &&
402     HcalIso<HcalIso_E &&
403 jueugste 1.2 missingHits<=missingHits_E) passed=true;
404 jueugste 1.1 // ECAL isolation
405     else if(combination[5]==1 &&
406     See<See_E &&
407     Dphi<Dphi_E &&
408     Deta<Deta_E &&
409     HoE<HoE_E &&
410     TrkIso<TrkIso_E &&
411     HcalIso<HcalIso_E &&
412 jueugste 1.2 missingHits<=missingHits_E) passed=true;
413 jueugste 1.1 // HCAL isolation
414     else if(combination[6]==1 &&
415     See<See_E &&
416     Dphi<Dphi_E &&
417     Deta<Deta_E &&
418     HoE<HoE_E &&
419     TrkIso<TrkIso_E &&
420     EcalIso<EcalIso_E &&
421 jueugste 1.2 missingHits<=missingHits_E) passed=true;
422 jueugste 1.1 // Conversions
423     else if(combination[7]==1 &&
424     See<See_E &&
425     Dphi<Dphi_E &&
426     Deta<Deta_E &&
427     HoE<HoE_E &&
428     TrkIso<TrkIso_E &&
429     EcalIso<EcalIso_E &&
430     HcalIso<HcalIso_E) passed=true;
431    
432    
433    
434    
435    
436     } // endcaps
437     return passed;
438     }
439    
440 jueugste 1.2
441     testEvent baseClass::fillTestEvent() {
442     // fill the event struct
443     testEvent event;
444     for(int i=0; i<NEles;i++) {
445     TLorentzVector vec;
446     vec.SetPxPyPzE(ElPx[i],ElPy[i],ElPz[i],ElE[i]);
447     event.fourMomentum.push_back(vec);
448     }
449     return event;
450     }
451    
452     int baseClass::ElectronPreselection(int i) {
453     // function to make electron preselection
454     // 0: preselected electron
455     // 1: not accepted in eta
456     // 2: not accepted in pt
457     // 3: not accepted in H/E
458     // 4: not accepted in dphi
459     // 5: not accepted in deta
460    
461     //double etaMax_B=1.44;
462     double etaMax_E=2.5;
463 jueugste 1.3 double ptMin=10;
464 jueugste 1.2 double HoverEMin=0.15;
465     double deltaPhiMax=0.15;
466     double deltaEtaMax=0.02;
467    
468     bool cutOn_Eta = true;
469     bool cutOn_Pt = true;
470     bool cutOn_HoE = true;
471     bool cutOn_dPhi= true;
472     bool cutOn_dEta= true;
473    
474     double eta = ElSCEta[i];
475     // double eta = ElEta[i];
476    
477     // double dPhiCorr=dphiCorrections(ElEta[i], ElPhi[i]);
478     // double dEtaCorr=detaCorrections(ElEta[i], ElPhi[i]);
479     double dPhiCorr=dphiCorrections(eta, ElPhi[i]);
480     double dEtaCorr=detaCorrections(eta, ElPhi[i]);
481    
482     TVector3 track;
483     track.SetXYZ(ElPx[i], ElPy[i], ElPz[i]);
484     double trackEta=track.Eta();
485     // double pt=ElCaloEnergy[i]/cosh(trackEta);
486     double pt = ElPt[i];
487    
488     double hoe = ElHcalOverEcal[i];
489     double dphi = ElDeltaPhiSuperClusterAtVtx[i] - dPhiCorr;
490     double deta = ElDeltaEtaSuperClusterAtVtx[i] - dEtaCorr;
491     // double dphi = ElDeltaPhiSuperClusterAtVtx[i];
492     // double deta = ElDeltaEtaSuperClusterAtVtx[i];
493    
494     if(cutOn_Eta) if(fabs(eta) > etaMax_E) return 1;
495     if(cutOn_Pt) if(pt < ptMin) return 2;
496     if(cutOn_HoE) if(hoe > HoverEMin) return 3;
497     if(cutOn_dPhi) if(fabs(dphi) > deltaPhiMax) return 4;
498     if(cutOn_dEta) if(fabs(deta) > deltaEtaMax) return 5;
499    
500     return 0;
501     }
502    
503 jueugste 1.3 int baseClass::ETHElectronID(int i, NminusOneCutLabel c) {
504     // ETH ID
505     // returns integer values corresponding
506     // to the passed ID cuts
507     // 0 fully IDed and isolated
508     // 1 delta phi cut not passed
509     // 2 delta eta cut not passed
510     // 3 |1/E-1/p| cut not passed
511     // 4 isolation not passed
512     // 5 is electron from conversion
513    
514     double etaMax_B = 1.4442;
515     double etaMin_E = 1.5660;
516     double etaMax_E = 2.5000;
517     double dPhiMax_B = 0.02;
518     double dPhiMax_E = 0.02;
519     double dEtaMax_B = 0.004;
520     double dEtaMax_E = 0.006;
521     double oEoPMax_B = 0.005;
522     double oEoPMax_E = 0.007;
523     double isoMax_B = 0.1;
524     double isoMax_E = 0.1;
525    
526     bool cutOn_dPhi = true;
527     bool cutOn_dEta = true;
528     bool cutOn_oEoP = true;
529     bool cutOn_iso = true;
530     bool cutOn_conversion = true;
531    
532     double eta = ElSCEta[i];
533     //double pt = ElPt[i];
534     double dPhiCorr = dphiCorrections(eta, ElPhi[i]);
535     double dEtaCorr = detaCorrections(eta, ElPhi[i]);
536     double dPhi = fabs(ElDeltaPhiSuperClusterAtVtx[i]-dPhiCorr);
537     double dEta = fabs(ElDeltaEtaSuperClusterAtVtx[i]-dEtaCorr);
538    
539     double p = ElTrkMomAtVtx[i];
540     double e = ElCaloEnergy[i];
541     double oEoP = fabs( 1/e - 1/p );
542    
543     double trkIso = ElDR03TkSumPt[i];
544     double scEt = e * sin(ElTheta[i]);
545     double iso = trkIso / scEt;
546    
547     int hits=1;
548     double dist=0.02;
549     double cot=0.02;
550     bool isConversion;
551     if((fabs(ElConvPartnerTrkDist[i])<dist && fabs(ElConvPartnerTrkDCot[i])<cot) /*|| ElNumberOfMissingInnerHits[i]>hits*/) isConversion=true;
552     else isConversion=false;
553    
554     if(c==dphi) cutOn_dPhi = false;
555     if(c==deta) cutOn_dEta = false;
556     if(c==oeop) cutOn_oEoP = false;
557     if(c==combiso) cutOn_iso = false;
558     if(c==conversion) cutOn_conversion = false;
559    
560     if(fabs(eta)<=etaMax_B) {
561     if(cutOn_dPhi) if(dPhi > dPhiMax_B) return 1;
562     if(cutOn_dEta) if(dEta > dEtaMax_B) return 2;
563     if(cutOn_oEoP) if(oEoP > oEoPMax_B) return 3;
564     if(cutOn_iso) if(iso > isoMax_B) return 4;
565     if(cutOn_conversion) if(isConversion==true) return 5;
566     return 0;
567     }
568     if(fabs(eta)>etaMin_E && fabs(eta)<etaMax_E) {
569     if(cutOn_dPhi) if(dPhi > dPhiMax_E) return 1;
570     if(cutOn_dEta) if(dEta > dEtaMax_E) return 2;
571     if(cutOn_oEoP) if(oEoP > oEoPMax_E) return 3;
572     if(cutOn_iso) if(iso > isoMax_E) return 4;
573     if(cutOn_conversion) if(isConversion==true) return 5;
574     return 0;
575     }
576     else {
577     cout<<"WPelectronID: electron not in ECAL acceptance region!!"<<endl;
578     return -1;
579     }
580    
581     }
582    
583 jueugste 1.2 int baseClass::WPElectronID(int i, NminusOneCutLabel c, vector<double> cuts) {
584     // function to ID and isolate electrons
585     // return values (int):
586     // 0 totally IDed & isolated electron
587     // 1 See cut not passed
588     // 2 dEta cut not passed
589     // 3 dPhi
590     // 4 HoE
591     // 5 combIso
592     // 6 conversion
593    
594 jueugste 1.3 double etaMax_B=1.4442;
595     double etaMin_E=1.5660;
596     double etaMax_E=2.5000;
597 jueugste 1.2 double SeeMax_B=cuts[0];
598     double dPhiMax_B=cuts[1];
599     double dEtaMax_B=cuts[2];
600     double HoEMax_B=cuts[3];
601     double TrkIsoMax_B=cuts[4];
602     double EcalIsoMax_B=cuts[5];
603     double HcalIsoMax_B=cuts[6];
604     double combIsoMax_B=cuts[7];
605     double SeeMax_E=cuts[8];
606     double dPhiMax_E=cuts[9];
607     double dEtaMax_E=cuts[10];
608     double HoEMax_E=cuts[11];
609     double TrkIsoMax_E=cuts[12];
610     double EcalIsoMax_E=cuts[13];
611     double HcalIsoMax_E=cuts[14];
612     double combIsoMax_E=cuts[15];
613     int hits=cuts[16];
614     double dist=cuts[17];
615     double cot=dist;
616    
617     bool cutOn_See=true;
618     bool cutOn_dPhi=true;
619     bool cutOn_dEta=true;
620     bool cutOn_HoE=true;
621     bool cutOn_combIso=true;
622     bool cutOn_conversion=true;
623    
624     // check geometrical acceptance
625     // for the isolation use the SC eta instead
626     // of the electron eta
627     double eta = ElSCEta[i];
628     // double eta = ElEta[i];
629    
630 jueugste 1.6 // double dPhiCorr=dphiCorrections(eta, ElPhi[i]);
631     // double dEtaCorr=detaCorrections(eta, ElPhi[i]);
632     double dPhiCorr=0.;//dphiCorrections(eta, ElPhi[i]);
633     double dEtaCorr=0.;//detaCorrections(eta, ElPhi[i]);
634 jueugste 1.2
635     double See = ElSigmaIetaIeta[i];
636     // apply the correction for the endcaps
637     double dPhi = fabs(ElDeltaPhiSuperClusterAtVtx[i]-dPhiCorr);
638     double dEta = fabs(ElDeltaEtaSuperClusterAtVtx[i]-dEtaCorr);
639     double HoE = ElHcalOverEcal[i];
640     double trkIso = ElDR03TkSumPt[i];
641     double ecalIso = ElDR03EcalRecHitSumEt[i];
642     double hcalIso = ElDR03HcalTowerSumEt[i];
643     double combIso_B = (trkIso + max(0., ecalIso - 1.) + hcalIso) / ElPt[i];
644     double combIso_E = (trkIso + ecalIso + hcalIso) / ElPt[i];
645     bool isConversion;
646 jueugste 1.6 if((fabs(ElConvPartnerTrkDist[i])<=dist && fabs(ElConvPartnerTrkDCot[i])<=cot) || ElNumberOfMissingInnerHits[i]>hits) isConversion=true;
647 jueugste 1.2 else isConversion=false;
648    
649     if(c==see) cutOn_See = false;
650     if(c==dphi) cutOn_dPhi = false;
651     if(c==deta) cutOn_dEta = false;
652     if(c==hoe) cutOn_HoE = false;
653     if(c==combiso) cutOn_combIso = false;
654     if(c==conversion) cutOn_conversion = false;
655    
656     if(fabs(eta)<=etaMax_B) {
657 jueugste 1.6 if(ElIDsimpleWP80relIso[i] < 5 || ElIDsimpleWP80relIso[i] == 6) return 1;
658     // if(cutOn_See) if(See > SeeMax_B) return 1;
659     // if(cutOn_dPhi) if(dPhi > dPhiMax_B) return 2;
660     // if(cutOn_dEta) if(dEta > dEtaMax_B) return 3;
661     // if(cutOn_HoE) if(HoE > HoEMax_B) return 4;
662 jueugste 1.2 if(cutOn_combIso) if(combIso_B > combIsoMax_B) return 5;
663 jueugste 1.6 // if(cutOn_combIso) if(trkIso / ElPt[i] < 0.09 && ecalIso / ElPt[i] < 0.07 && hcalIso / ElPt[i] < 0.1) return 5;
664     // if(cutOn_conversion) if(isConversion==true) return 6;
665 jueugste 1.2 return 0;
666     }
667     if(fabs(eta)>etaMin_E && fabs(eta)<etaMax_E) {
668 jueugste 1.6 if(ElIDsimpleWP80relIso[i] < 5 || ElIDsimpleWP80relIso[i] == 6) return 1;
669     // if(cutOn_See) if(See > SeeMax_E) return 1;
670     // if(cutOn_dPhi) if(dPhi > dPhiMax_E) return 2;
671     // if(cutOn_dEta) if(dEta > dEtaMax_E) return 3;
672     // if(cutOn_HoE) if(HoE > HoEMax_E) return 4;
673 jueugste 1.2 if(cutOn_combIso) if(combIso_E > combIsoMax_E) return 5;
674 jueugste 1.6 // if(cutOn_combIso) if(trkIso / ElPt[i] < 0.04 && ecalIso / ElPt[i] < 0.05 && hcalIso / ElPt[i] < 0.025) return 5;
675     // if(cutOn_conversion) if(isConversion==true) return 6;
676 jueugste 1.2 return 0;
677     }
678     else {
679 jueugste 1.3 // cout<<"WPelectronID: electron not in ECAL acceptance region!!"<<endl;
680 jueugste 1.2 return -1;
681     }
682     }
683    
684     double baseClass::detaCorrections(double etaEle,double phiEle) {
685     // corrects for misalignement in the endcpas
686     TF2 f12_fit2("f12_fit2","[0]+(TMath::TanH(y)/325.)*([1]-([2]*TMath::SinH(y)*TMath::Cos(x-[3])))",-TMath::Pi(),TMath::Pi(),-10.,10.);
687     if (etaEle>1.479)
688     {
689     f12_fit2.SetParameter(0,0.0013);
690     f12_fit2.SetParameter(1,-0.06);
691     f12_fit2.SetParameter(2,0.52);
692     f12_fit2.SetParameter(3,2.17);
693     }
694     else if (etaEle<-1.479)
695     {
696     f12_fit2.SetParameter(0,-0.0013);
697     f12_fit2.SetParameter(1,-0.32);
698     f12_fit2.SetParameter(2,0.45);
699     f12_fit2.SetParameter(3,-1.58);
700     }
701     return f12_fit2.Eval(phiEle,etaEle);
702     }
703    
704     double baseClass::dphiCorrections(double etaEle,double phiEle) {
705     // corrects for misalignement in the endcpas
706     TF2 f12_fit2("f12_fit2","[0]+[1]*(TMath::SinH(y)/325.)*(TMath::Sin([2]-x))",-TMath::Pi(),TMath::Pi(),-10.,10.);
707     if (etaEle>1.479)
708     {
709     f12_fit2.SetParameter(0,0.0);
710     f12_fit2.SetParameter(1,0.52);
711     f12_fit2.SetParameter(2,2.17);
712     }
713     else if (etaEle<-1.479)
714     {
715     f12_fit2.SetParameter(0,0.0);
716     f12_fit2.SetParameter(1,0.45);
717     f12_fit2.SetParameter(2,-1.58);
718     }
719     return f12_fit2.Eval(phiEle,etaEle);
720     }
721    
722    
723     int baseClass::WSelection(int i, WNminusOneCutLabel c) {
724     // idea: write cuts into a vector and give the vector to the function
725     // double etaMax_B=1.44;
726     double pre_metMin=20;
727     double pre_ptMin=20;
728     double metMin=30;
729     double ptMin=25;
730     double mtMin=60;
731    
732     bool cutOn_Preselection=true;
733     bool cutOn_Met=true;
734     bool cutOn_Pt=true;
735     bool cutOn_Mt=true;
736    
737     double met=PFMET;
738     double pt=ElPt[i];
739     double delPhiMet=deltaPhi(ElPhi[i], PFMETphi);
740     double mt=transverseMass(pt, met, delPhiMet);
741     // check geometrical acceptance
742     // for the isolation use the SC eta instead
743     // of the electron eta
744     // double eta = ElSCEta[i];
745     // double eta = ElEta[i];
746    
747     if(c==pre) cutOn_Preselection = false;
748     if(c==misset) cutOn_Met = false;
749     if(c==pt) cutOn_Pt = false;
750     if(c==mt) cutOn_Mt = false;
751    
752    
753     // if(fabs(eta)<=etaMax_B) {
754     if(cutOn_Preselection) if(met<pre_metMin || pt<pre_ptMin) return 1;
755     if(cutOn_Met) if(met<metMin) return 2;
756     if(cutOn_Pt) if(pt<ptMin) return 3;
757     if(cutOn_Mt) if(mt<mtMin) return 4;
758     return 0;
759     // }
760     // if(fabs(eta)>etaMax_B) {
761     // if(cutOn_Preselection) if(met<pre_metMin || pt<pre_ptMin) return 1;
762     // if(cutOn_Met) if(met<metMin) return 2;
763     // if(cutOn_Pt) if(pt<ptMin) return 3;
764     // if(cutOn_Mt) if(mt<mtMin) return 4;
765     // return 0;
766     // }
767     }
768    
769    
770    
771     // void initializeWP(vector<double>cuts) {
772     // SeeMax_B=cuts[0];
773     // DphiMax_B=cuts[1];
774     // DetaMax_B=cuts[2];
775     // HoEMax_B=cuts[3];
776     // combIsoMax_B=cuts[4];
777     // // TrkIsoMax_B=cuts[4];
778     // // EcalIsoMax_B=cuts[5];
779     // // HcalIsoMax_B=cuts[6];
780     // SeeMax_E=cuts[8];
781     // DphiMax_E=cuts[9];
782     // DetaMax_E=cuts[10];
783     // HoEMax_E=cuts[11];
784     // // TrkIsoMax_E=cuts[12];
785     // // EcalIsoMax_E=cuts[13];
786     // // HcalIsoMax_E=cuts[14];
787     // combIsoMax_E=cuts[12];
788    
789     // }
790    
791 jueugste 1.1 bool baseClass::ethNminus1(int isECALdriven, double deltaPhi, double deltaEta, double e, double p, double eta, char *cut, bool iso)
792     {
793     bool passed=false;
794     double oneOverEoneOverP=1/e-1/p;
795     // cut values
796     double etaMax_B=1.44;
797     double etaMax_E=2.6;
798     double deltaPhiCut_B=0.02;
799     double deltaPhiCut_E=0.02;
800     double deltaEtaCut_B=0.004;
801     double deltaEtaCut_E=0.006;
802     double oneOverEoneOverPCut_B=0.005;
803     double oneOverEoneOverPCut_E=0.007;
804     int combination [4] = {0,0,0,0};
805     if(!strcmp(cut,"deltaPhi")) combination [0] = 1;
806     else if(!strcmp(cut,"deltaEta")) combination [1] =1;
807     else if(!strcmp(cut,"oneOverEminusOneOverP")) combination [2] = 1;
808     else if(!strcmp(cut,"iso")) combination[3] = 1;
809     else {
810     cout<<"ERROR: not able to give N-1 flag for the ETH!"<<endl;
811     }
812     if(isECALdriven == 1) {
813     // barrel
814     if(fabs(eta)<etaMax_B) {
815     // delta phi
816     if(combination[0]==1 && fabs(deltaEta)<deltaEtaCut_B && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_B && iso==true) {
817     passed = true;
818     }
819     // delta eta
820     else if(combination[1]==1 && fabs(deltaPhi)<deltaPhiCut_B && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_B && iso==true) {
821     passed = true;
822     }
823    
824    
825     // |1/E-1/p|
826     else if(combination[2]==1 && fabs(deltaPhi)<deltaPhiCut_B && fabs(deltaEta)< deltaEtaCut_B && iso==true) {
827     passed = true;
828     }
829     // iso
830     else if((combination[3]==1 && fabs(deltaPhi)<deltaPhiCut_B && fabs(deltaEta)< deltaEtaCut_B && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_B)) {
831     passed=true;
832     }
833     }
834     //endcap
835     else if(fabs(eta)>etaMax_B && fabs(eta)<etaMax_E) {
836     // delta phi
837     if(combination[0]==1 && fabs(deltaEta)< deltaEtaCut_E && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_E && iso==true) {
838     passed = true;
839     }
840     // delta eta
841     else if(combination[1]==1 && fabs(deltaPhi)< deltaPhiCut_E && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_E && iso==true) {
842     passed = true;
843     }
844     // |1/E-1/p|
845     else if(combination[2]==1 && fabs(deltaPhi)< deltaPhiCut_E && fabs(deltaEta)< deltaEtaCut_E && iso ==true) {
846     passed = true;
847     }
848     // iso
849     else if((combination[3]==1 && fabs(deltaPhi)< deltaPhiCut_E && fabs(deltaEta)< deltaEtaCut_E && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_E)) {
850     passed=true;
851     }
852     } // endcap
853     } // ecal driven
854     return passed;
855     }
856    
857    
858     bool baseClass::ethID(int isECALdriven, double deltaPhi, double deltaEta, double e, double p, double eta)
859     {
860     bool passed=false;
861     double oneOverEoneOverP=1/e-1/p;
862     // cut values
863     double etaMax_B=1.44;
864     double etaMax_E=2.6;
865     double deltaPhiCut_B=0.02;
866     double deltaPhiCut_E=0.02;
867     double deltaEtaCut_B=0.004;
868     double deltaEtaCut_E=0.006;
869     double oneOverEoneOverPCut_B=0.005;
870     double oneOverEoneOverPCut_E=0.007;
871     if(isECALdriven == 1 &&
872     fabs(eta)<etaMax_B &&
873     fabs(deltaPhi)< deltaPhiCut_B &&
874     fabs(deltaEta)< deltaEtaCut_B &&
875     fabs(oneOverEoneOverP)<oneOverEoneOverPCut_B) {
876     passed=true;
877     }
878     else if(isECALdriven == 1 &&
879     fabs(eta)>etaMax_B &&
880     fabs(eta)<etaMax_E &&
881     fabs(deltaPhi)< deltaPhiCut_E &&
882     fabs(deltaEta)< deltaEtaCut_E &&
883     fabs(oneOverEoneOverP)<oneOverEoneOverPCut_E) {
884     passed=true;
885     }
886     return passed;
887     }
888    
889     bool baseClass::isolation(double sumpt, double e)
890     {
891     bool passed =false;
892     double iso;
893     iso = sumpt/e;
894     if(iso<0.1) {
895     passed =true;
896     }
897     return passed;
898     }
899    
900     double baseClass::deltaPhi(double phi1, double phi2) {
901     double delPhi;
902     double Pi=2*acos(0);
903     delPhi=fabs(phi1-phi2);
904     if(delPhi>Pi) {
905     delPhi=2*Pi-delPhi;
906     }
907     return delPhi;
908     }
909    
910     double baseClass::transverseMass(double pt1, double pt2, double delPhi) {
911     // pt2=met!!
912     double mt;
913     mt=sqrt(2*fabs(pt1)*fabs(pt2)*(1-cos(delPhi)));
914     return mt;
915     }
916    
917     void baseClass::init()
918     {
919     STDOUT("begin");
920     //std::cout << "baseClass::init(): begin " << std::endl;
921    
922     tree_ = NULL;
923     readInputList();
924 jueugste 1.2 // readCutFile(); // not using cut file, need dummy cut file as input...
925 jueugste 1.1 if(tree_ == NULL){
926     std::cout << "baseClass::init(): ERROR: tree_ = NULL " << std::endl;
927     return;
928     }
929     Init(tree_);
930    
931     //char output_root_title[200];
932     //sprintf(output_root_title,"%s%s",&std::string(*outputFileName_)[0],".root");
933     //output_root_ = new TFile(&output_root_title[0],"RECREATE");
934    
935     //directly from string
936     output_root_ = new TFile((*outputFileName_ + ".root").c_str(),"RECREATE");
937    
938     // for (map<string, cut>::iterator it = cutName_cut_.begin();
939     // it != cutName_cut_.end(); it++) STDOUT("cutName_cut->first = "<<it->first)
940     // for (vector<string>::iterator it = orderedCutNames_.begin();
941     // it != orderedCutNames_.end(); it++) STDOUT("orderedCutNames_ = "<<*it)
942     STDOUT("end");
943     }
944    
945     void baseClass::readInputList()
946     {
947    
948     TChain *chain = new TChain(treeName_->c_str());
949     char pName[500];
950     skimWasMade_ = true;
951     NBeforeSkim_ = 0;
952     int NBeforeSkim;
953    
954     std::cout << "baseClass::readinputList(): inputList_ = "<< *inputList_ << std::endl;
955    
956     ifstream is(inputList_->c_str());
957     if(is.good())
958     {
959     cout << "baseClass::readInputList: Reading list: " << *inputList_ << " ......." << endl;
960     while( is.getline(pName, 500, '\n') )
961     {
962     if (pName[0] == '#') continue;
963     //if (pName[0] == ' ') continue; // do we want to skip lines that start with a space?
964     if (pName[0] == '\n') continue;// simple protection against blank lines
965     STDOUT("Adding file: " << pName);
966     chain->Add(pName);
967     NBeforeSkim = getGlobalInfoNstart(pName);
968     NBeforeSkim_ = NBeforeSkim_ + NBeforeSkim;
969     STDOUT("Initial number of events: NBeforeSkim, NBeforeSkim_ = "<<NBeforeSkim<<", "<<NBeforeSkim_);
970     }
971 jueugste 1.2
972 jueugste 1.1 tree_ = chain;
973     cout << "baseClass::readInputList: Finished reading list: " << *inputList_ << endl;
974     }
975     else
976     {
977     cout << "baseClass::readInputList: ERROR opening inputList:" << *inputList_ << endl;
978     exit (1);
979     }
980     is.close();
981    
982     }
983    
984     void baseClass::readCutFile()
985     {
986     string s;
987     STDOUT("Reading cutFile_ = "<< *cutFile_)
988    
989     ifstream is(cutFile_->c_str());
990     if(is.good())
991     {
992     // STDOUT("Reading file: " << *cutFile_ );
993     int id=0;
994     int optimize_count=0;
995     while( getline(is,s) )
996     {
997     STDOUT("read line: " << s);
998     if (s[0] == '#' || s.empty()) continue;
999     vector<string> v = split(s);
1000     if ( v.size() == 0 ) continue;
1001     if (v[1]=="OPT") // add code for grabbing optimizer objects
1002     {
1003     if (optimizeName_cut_.size()>=6)
1004     {
1005     STDOUT("WARNING: Optimizer can only accept up to 6 variables.\nVariable "<<v[0]<<" is not being included.");
1006     continue;
1007     }
1008     bool found=false;
1009     for (int i=0;i<optimizeName_cut_.size();++i)
1010     {
1011     if (optimizeName_cut_[i].variableName==v[0])
1012     {
1013     STDOUT("ERROR: variableName = "<<v[0]<<" is already being optimized in optimizedName_cut_. Skipping.");
1014     found=true;
1015     break;
1016     }
1017     }
1018     if (found) continue;
1019    
1020     int level_int = atoi(v[5].c_str());
1021     bool greaterthan=true;
1022     if (v[2]=="<") greaterthan=false;
1023     double minval=atof(v[3].c_str());
1024     double maxval=atof(v[4].c_str());
1025     Optimize opt(optimize_count,v[0],minval, maxval, greaterthan, level_int);
1026     optimizeName_cut_[optimize_count]=opt; // order cuts by cut #, rather than name, so that optimization histogram is consistently ordered
1027     ++optimize_count;
1028     continue;
1029     }
1030    
1031     map<string, cut>::iterator cc = cutName_cut_.find(v[0]);
1032     if( cc != cutName_cut_.end() )
1033     {
1034     STDOUT("ERROR: variableName = "<< v[0] << " exists already in cutName_cut_. Returning.");
1035     return;
1036     }
1037    
1038     int level_int = atoi( v[5].c_str() );
1039     if(level_int == -1)
1040     {
1041     map<string, preCut>::iterator cc = preCutName_cut_.find(v[0]);
1042     if( cc != preCutName_cut_.end() )
1043     {
1044     STDOUT("ERROR: variableName = "<< v[0] << " exists already in preCutName_cut_. Returning.");
1045     return;
1046     }
1047     preCutInfo_ << "### Preliminary cut values: " << s <<endl;
1048     preCut thisPreCut;
1049     thisPreCut.variableName = v[0];
1050     thisPreCut.value1 = decodeCutValue( v[1] );
1051     thisPreCut.value2 = decodeCutValue( v[2] );
1052     thisPreCut.value3 = decodeCutValue( v[3] );
1053     thisPreCut.value4 = decodeCutValue( v[4] );
1054     preCutName_cut_[thisPreCut.variableName]=thisPreCut;
1055     continue;
1056     }
1057     cut thisCut;
1058     thisCut.variableName = v[0];
1059     string m1=v[1];
1060     string M1=v[2];
1061     string m2=v[3];
1062     string M2=v[4];
1063     if( m1=="-" || M1=="-" )
1064     {
1065     STDOUT("ERROR: minValue1 and maxValue2 have to be provided. Returning.");
1066     return; // FIXME implement exception
1067     }
1068     if( (m2=="-" && M2!="-") || (m2!="-" && M2=="-") )
1069     {
1070     STDOUT("ERROR: if any of minValue2 and maxValue2 is -, then both have to be -. Returning");
1071     return; // FIXME implement exception
1072     }
1073     if( m2=="-") m2="+inf";
1074     if( M2=="-") M2="-inf";
1075     thisCut.minValue1 = decodeCutValue( m1 );
1076     thisCut.maxValue1 = decodeCutValue( M1 );
1077     thisCut.minValue2 = decodeCutValue( m2 );
1078     thisCut.maxValue2 = decodeCutValue( M2 );
1079     thisCut.level_int = level_int;
1080     thisCut.level_str = v[5];
1081     thisCut.histoNBins = atoi( v[6].c_str() );
1082     thisCut.histoMin = atof( v[7].c_str() );
1083     thisCut.histoMax = atof( v[8].c_str() );
1084     // Not filled from file
1085     thisCut.id=++id;
1086     string s1;
1087     if(skimWasMade_)
1088     {
1089     s1 = "cutHisto_skim___________________" + thisCut.variableName;
1090     }
1091     else
1092     {
1093     s1 = "cutHisto_noCuts_________________" + thisCut.variableName;
1094     }
1095     string s2 = "cutHisto_allPreviousCuts________" + thisCut.variableName;
1096     string s3 = "cutHisto_allOthrSmAndLwrLvlCuts_" + thisCut.variableName;
1097     string s4 = "cutHisto_allOtherCuts___________" + thisCut.variableName;
1098     string s5 = "cutHisto_allCuts________________" + thisCut.variableName;
1099     thisCut.histo1 = TH1F (s1.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
1100     thisCut.histo2 = TH1F (s2.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
1101     thisCut.histo3 = TH1F (s3.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
1102     thisCut.histo4 = TH1F (s4.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
1103     thisCut.histo5 = TH1F (s5.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
1104     thisCut.histo1.Sumw2();
1105     thisCut.histo2.Sumw2();
1106     thisCut.histo3.Sumw2();
1107     thisCut.histo4.Sumw2();
1108     thisCut.histo5.Sumw2();
1109     // Filled event by event
1110     thisCut.filled = false;
1111     thisCut.value = 0;
1112     thisCut.passed = false;
1113     thisCut.nEvtInput=0;
1114     thisCut.nEvtPassed=0;
1115    
1116     orderedCutNames_.push_back(thisCut.variableName);
1117     cutName_cut_[thisCut.variableName]=thisCut;
1118    
1119     }
1120     STDOUT( "baseClass::readCutFile: Finished reading cutFile: " << *cutFile_ );
1121     }
1122     else
1123     {
1124     STDOUT("ERROR opening cutFile:" << *cutFile_ );
1125     exit (1);
1126     }
1127     // make optimizer histogram
1128     if (optimizeName_cut_.size()>0)
1129     {
1130     h_optimizer_=new TH1F("optimizer","Optimization of cut variables",(int)pow(10,optimizeName_cut_.size()),0,
1131     pow(10,optimizeName_cut_.size()));
1132     }
1133    
1134     // Create a histogram that will show events passing cuts
1135     int cutsize=orderedCutNames_.size()+1;
1136     if (skimWasMade_) ++cutsize;
1137     gDirectory->cd();
1138     eventcuts_=new TH1F("EventsPassingCuts","Events Passing Cuts",cutsize,0,cutsize);
1139    
1140     is.close();
1141    
1142     }
1143    
1144     void baseClass::resetCuts()
1145     {
1146     for (map<string, cut>::iterator cc = cutName_cut_.begin(); cc != cutName_cut_.end(); cc++)
1147     {
1148     cut * c = & (cc->second);
1149     c->filled = false;
1150     c->value = 0;
1151     c->passed = false;
1152     }
1153     combCutName_passed_.clear();
1154     return;
1155     }
1156    
1157     void baseClass::fillVariableWithValue(const string& s, const double& d)
1158     {
1159     map<string, cut>::iterator cc = cutName_cut_.find(s);
1160     if( cc == cutName_cut_.end() )
1161     {
1162     //STDOUT("variableName = "<< s << " not found in cutName_cut_.");
1163     return;
1164     }
1165     else
1166     {
1167     cut * c = & (cc->second);
1168     c->filled = true;
1169     c->value = d;
1170     }
1171     fillOptimizerWithValue(s, d);
1172     return;
1173     }
1174    
1175     bool baseClass::variableIsFilled(const string& s)
1176     {
1177     map<string, cut>::iterator cc = cutName_cut_.find(s);
1178     if( cc == cutName_cut_.end() )
1179     {
1180     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
1181     }
1182     cut * c = & (cc->second);
1183     return (c->filled);
1184     }
1185    
1186     double baseClass::getVariableValue(const string& s)
1187     {
1188     map<string, cut>::iterator cc = cutName_cut_.find(s);
1189     if( cc == cutName_cut_.end() )
1190     {
1191     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_.");
1192     }
1193     cut * c = & (cc->second);
1194     if( !variableIsFilled(s) )
1195     {
1196     STDOUT("ERROR: requesting value of not filled variable "<<s);
1197     }
1198     return (c->value);
1199     }
1200    
1201     void baseClass::fillOptimizerWithValue(const string& s, const double& d)
1202     {
1203     for (int i=0;i<optimizeName_cut_.size();++i)
1204     {
1205     if (optimizeName_cut_[i].variableName==s)
1206     {
1207     optimizeName_cut_[i].value=d;
1208     return;
1209     }
1210     }
1211     return;
1212     }
1213    
1214     void baseClass::evaluateCuts()
1215     {
1216     // resetCuts();
1217     combCutName_passed_.clear();
1218     for (vector<string>::iterator it = orderedCutNames_.begin();
1219     it != orderedCutNames_.end(); it++)
1220     {
1221     cut * c = & (cutName_cut_.find(*it)->second);
1222     if( ! ( c->filled && (c->minValue1 < c->value && c->value <= c->maxValue1 || c->minValue2 < c->value && c->value <= c->maxValue2 ) ) )
1223     {
1224     c->passed = false;
1225     combCutName_passed_[c->level_str.c_str()] = false;
1226     combCutName_passed_["all"] = false;
1227     }
1228     else
1229     {
1230     c->passed = true;
1231     map<string,bool>::iterator cp = combCutName_passed_.find( c->level_str.c_str() );
1232     combCutName_passed_[c->level_str.c_str()] = (cp==combCutName_passed_.end()?true:cp->second);
1233     map<string,bool>::iterator ap = combCutName_passed_.find( "all" );
1234     combCutName_passed_["all"] = (ap==combCutName_passed_.end()?true:ap->second);
1235     }
1236     }
1237    
1238     // reset optimization cut values
1239     //for (int i=0;i<optimizeName_cut_.size();++i)
1240     // optimizeName_cut_[i].value=0;
1241     runOptimizer();
1242    
1243     if( !fillCutHistos() )
1244     {
1245     STDOUT("ERROR: fillCutHistos did not complete successfully.");
1246     }
1247    
1248     if( !updateCutEffic() )
1249     {
1250     STDOUT("ERROR: updateCutEffic did not complete successfully.");
1251     }
1252    
1253     return ;
1254     }
1255    
1256     void baseClass::runOptimizer()
1257     {
1258    
1259     // don't run optimizer if no optimized cuts specified
1260     if (optimizeName_cut_.size()==0)
1261     return;
1262    
1263     // first, check that all cuts (except those to be optimized) have been passed
1264    
1265     for (vector<string>::iterator it = orderedCutNames_.begin();
1266     it != orderedCutNames_.end(); it++)
1267     {
1268     bool ignorecut=false;
1269     for (unsigned int i=0; i < optimizeName_cut_.size();++i)
1270     {
1271     const string str = (const string)(*it);
1272     if (optimizeName_cut_[i].variableName.compare(str)==0)
1273     {
1274     ignorecut=true;
1275     break;
1276     }
1277     }
1278     if (ignorecut) continue;
1279     if (passedCut(*it) == false)
1280     return;
1281     }
1282    
1283     /*
1284     if (combCutName_passed_["all"] == false)
1285     return;
1286     */
1287    
1288     // loop over up to 6 cuts
1289     int counter=0;
1290     int thesize=optimizeName_cut_.size();
1291     int mysize=thesize;
1292     std::vector<bool> counterbins;
1293     for (int i=0;i<pow(10,thesize);++i) counterbins.push_back(true); // assume true
1294    
1295     // lowest-numbered cut appears first in cut ordering
1296     // That is, for cut: ABCDEF
1297     // A is the index of cut0, B is cut 1, etc.
1298     for (int cc=0;cc<thesize;++cc) // loop over all cuts, starting at cut 0
1299     {
1300     --mysize;
1301     for (int i=0;i<10;++i) // loop over 10 cuts for each
1302     {
1303     if (!optimizeName_cut_[cc].Compare(i)) // cut failed; set all values associated with cut to false
1304     {
1305     // loop over all cut values starting with current cut
1306     for (unsigned int j=(int)(i*pow(10,mysize));j<int(pow(10,thesize));++j)
1307     {
1308     // if relevant digit of the cut value matches the current (failed) cut, set this cut to false
1309     if ((j/int(pow(10,mysize)))%10==i)
1310     counterbins[j]=false;
1311     if (j>counterbins.size())
1312     continue; // shouldn't ever happen
1313     }
1314     } // if (cut comparison failed)
1315     } // for (int i=0;i<10;++i)
1316     }
1317     // now fill histograms
1318     for (int i=0;i<counterbins.size();++i)
1319     {
1320     if (counterbins[i]==true)
1321     h_optimizer_->Fill(i,1);
1322    
1323     }
1324    
1325     return;
1326     } //runOptimizer
1327    
1328     bool baseClass::passedCut(const string& s)
1329     {
1330     bool ret = false;
1331     // map<string, bool>::iterator vp = cutName_passed_.find(s);
1332     // if( vp != cutName_passed_.end() ) return ret = vp->second;
1333     map<string, cut>::iterator cc = cutName_cut_.find(s);
1334     if( cc != cutName_cut_.end() )
1335     {
1336     cut * c = & (cutName_cut_.find(s)->second);
1337     return (c->filled && c->passed);
1338     }
1339     map<string, bool>::iterator cp = combCutName_passed_.find(s);
1340     if( cp != combCutName_passed_.end() )
1341     {
1342     return ret = cp->second;
1343     }
1344     STDOUT("ERROR: did not find variableName = "<<s<<" neither in cutName_cut_ nor combCutName_passed_. Returning false.");
1345     return (ret=false);
1346     }
1347    
1348     bool baseClass::passedAllPreviousCuts(const string& s)
1349     {
1350     //STDOUT("Examining variableName = "<<s);
1351    
1352     map<string, cut>::iterator cc = cutName_cut_.find(s);
1353     if( cc == cutName_cut_.end() )
1354     {
1355     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1356     return false;
1357     }
1358    
1359     for (vector<string>::iterator it = orderedCutNames_.begin();
1360     it != orderedCutNames_.end(); it++)
1361     {
1362     cut * c = & (cutName_cut_.find(*it)->second);
1363     if( c->variableName == s )
1364     {
1365     return true;
1366     }
1367     else
1368     {
1369     if( ! (c->filled && c->passed) ) return false;
1370     }
1371     }
1372     STDOUT("ERROR. It should never pass from here.");
1373     }
1374    
1375     bool baseClass::passedAllOtherCuts(const string& s)
1376     {
1377     //STDOUT("Examining variableName = "<<s);
1378     bool ret = true;
1379    
1380     map<string, cut>::iterator cc = cutName_cut_.find(s);
1381     if( cc == cutName_cut_.end() )
1382     {
1383     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1384     return false;
1385     }
1386    
1387     for (map<string, cut>::iterator cc = cutName_cut_.begin(); cc != cutName_cut_.end(); cc++)
1388     {
1389     cut * c = & (cc->second);
1390     if( c->variableName == s )
1391     {
1392     continue;
1393     }
1394     else
1395     {
1396     if( ! (c->filled && c->passed) ) return false;
1397     }
1398     }
1399     return ret;
1400     }
1401    
1402     bool baseClass::passedAllOtherSameAndLowerLevelCuts(const string& s)
1403     {
1404     //STDOUT("Examining variableName = "<<s);
1405     bool ret = true;
1406     int cutLevel;
1407     map<string, cut>::iterator cc = cutName_cut_.find(s);
1408     if( cc == cutName_cut_.end() )
1409     {
1410     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1411     return false;
1412     }
1413     else
1414     {
1415     cutLevel = cc->second.level_int;
1416     }
1417    
1418     for (map<string, cut>::iterator cc = cutName_cut_.begin(); cc != cutName_cut_.end(); cc++)
1419     {
1420     cut * c = & (cc->second);
1421     if( c->level_int > cutLevel || c->variableName == s )
1422     {
1423     continue;
1424     }
1425     else
1426     {
1427     if( ! (c->filled && c->passed) ) return false;
1428     }
1429     }
1430     return ret;
1431     }
1432    
1433     double baseClass::getPreCutValue1(const string& s)
1434     {
1435     double ret;
1436     map<string, preCut>::iterator cc = preCutName_cut_.find(s);
1437     if( cc == preCutName_cut_.end() )
1438     {
1439     STDOUT("ERROR: did not find variableName = "<<s<<" in preCutName_cut_. Returning");
1440     }
1441     preCut * c = & (cc->second);
1442     if(c->value1 == -9999999) STDOUT("ERROR: value1 of preliminary cut "<<s<<" was not provided.");
1443     return (c->value1);
1444     }
1445    
1446     double baseClass::getPreCutValue2(const string& s)
1447     {
1448     double ret;
1449     map<string, preCut>::iterator cc = preCutName_cut_.find(s);
1450     if( cc == preCutName_cut_.end() )
1451     {
1452     STDOUT("ERROR: did not find variableName = "<<s<<" in preCutName_cut_. Returning");
1453     }
1454     preCut * c = & (cc->second);
1455     if(c->value2 == -9999999) STDOUT("ERROR: value2 of preliminary cut "<<s<<" was not provided.");
1456     return (c->value2);
1457     }
1458    
1459     double baseClass::getPreCutValue3(const string& s)
1460     {
1461     double ret;
1462     map<string, preCut>::iterator cc = preCutName_cut_.find(s);
1463     if( cc == preCutName_cut_.end() )
1464     {
1465     STDOUT("ERROR: did not find variableName = "<<s<<" in preCutName_cut_. Returning");
1466     }
1467     preCut * c = & (cc->second);
1468     if(c->value3 == -9999999) STDOUT("ERROR: value3 of preliminary cut "<<s<<" was not provided.");
1469     return (c->value3);
1470     }
1471    
1472     double baseClass::getPreCutValue4(const string& s)
1473     {
1474     double ret;
1475     map<string, preCut>::iterator cc = preCutName_cut_.find(s);
1476     if( cc == preCutName_cut_.end() )
1477     {
1478     STDOUT("ERROR: did not find variableName = "<<s<<" in preCutName_cut_. Returning");
1479     }
1480     preCut * c = & (cc->second);
1481     if(c->value4 == -9999999) STDOUT("ERROR: value4 of preliminary cut "<<s<<" was not provided.");
1482     return (c->value4);
1483     }
1484    
1485    
1486     double baseClass::getCutMinValue1(const string& s)
1487     {
1488     double ret;
1489     map<string, cut>::iterator cc = cutName_cut_.find(s);
1490     if( cc == cutName_cut_.end() )
1491     {
1492     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
1493     }
1494     cut * c = & (cc->second);
1495     return (c->minValue1);
1496     }
1497    
1498     double baseClass::getCutMaxValue1(const string& s)
1499     {
1500     double ret;
1501     map<string, cut>::iterator cc = cutName_cut_.find(s);
1502     if( cc == cutName_cut_.end() )
1503     {
1504     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
1505     }
1506     cut * c = & (cc->second);
1507     return (c->maxValue1);
1508     }
1509    
1510     double baseClass::getCutMinValue2(const string& s)
1511     {
1512     double ret;
1513     map<string, cut>::iterator cc = cutName_cut_.find(s);
1514     if( cc == cutName_cut_.end() )
1515     {
1516     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
1517     }
1518     cut * c = & (cc->second);
1519     return (c->minValue2);
1520     }
1521    
1522     double baseClass::getCutMaxValue2(const string& s)
1523     {
1524     double ret;
1525     map<string, cut>::iterator cc = cutName_cut_.find(s);
1526     if( cc == cutName_cut_.end() )
1527     {
1528     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
1529     }
1530     cut * c = & (cc->second);
1531     return (c->maxValue2);
1532     }
1533    
1534     const TH1F& baseClass::getHisto_noCuts_or_skim(const string& s)
1535     {
1536     map<string, cut>::iterator cc = cutName_cut_.find(s);
1537     if( cc == cutName_cut_.end() )
1538     {
1539     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1540     }
1541    
1542     cut * c = & (cutName_cut_.find(s)->second);
1543     return (c->histo1);
1544     }
1545    
1546     const TH1F& baseClass::getHisto_allPreviousCuts(const string& s)
1547     {
1548     map<string, cut>::iterator cc = cutName_cut_.find(s);
1549     if( cc == cutName_cut_.end() )
1550     {
1551     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1552     }
1553    
1554     cut * c = & (cutName_cut_.find(s)->second);
1555     return (c->histo2);
1556     }
1557    
1558     const TH1F& baseClass::getHisto_allOthrSmAndLwrLvlCuts(const string& s)
1559     {
1560     map<string, cut>::iterator cc = cutName_cut_.find(s);
1561     if( cc == cutName_cut_.end() )
1562     {
1563     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1564     }
1565    
1566     cut * c = & (cutName_cut_.find(s)->second);
1567     return (c->histo3);
1568     }
1569    
1570     const TH1F& baseClass::getHisto_allOtherCuts(const string& s)
1571     {
1572     map<string, cut>::iterator cc = cutName_cut_.find(s);
1573     if( cc == cutName_cut_.end() )
1574     {
1575     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1576     }
1577    
1578     cut * c = & (cutName_cut_.find(s)->second);
1579     return (c->histo4);
1580     }
1581    
1582     const TH1F& baseClass::getHisto_allCuts(const string& s)
1583     {
1584     map<string, cut>::iterator cc = cutName_cut_.find(s);
1585     if( cc == cutName_cut_.end() )
1586     {
1587     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1588     }
1589    
1590     cut * c = & (cutName_cut_.find(s)->second);
1591     return (c->histo5);
1592     }
1593    
1594     int baseClass::getHistoNBins(const string& s)
1595     {
1596     map<string, cut>::iterator cc = cutName_cut_.find(s);
1597     if( cc == cutName_cut_.end() )
1598     {
1599     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1600     }
1601     cut * c = & (cutName_cut_.find(s)->second);
1602     return (c->histoNBins);
1603     }
1604    
1605     double baseClass::getHistoMin(const string& s)
1606     {
1607     map<string, cut>::iterator cc = cutName_cut_.find(s);
1608     if( cc == cutName_cut_.end() )
1609     {
1610     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1611     }
1612     cut * c = & (cutName_cut_.find(s)->second);
1613     return (c->histoMin);
1614     }
1615    
1616     double baseClass::getHistoMax(const string& s)
1617     {
1618     map<string, cut>::iterator cc = cutName_cut_.find(s);
1619     if( cc == cutName_cut_.end() )
1620     {
1621     STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
1622     }
1623     cut * c = & (cutName_cut_.find(s)->second);
1624     return (c->histoMax);
1625     }
1626    
1627    
1628     bool baseClass::fillCutHistos()
1629     {
1630     bool ret = true;
1631     for (vector<string>::iterator it = orderedCutNames_.begin();
1632     it != orderedCutNames_.end(); it++)
1633     {
1634     cut * c = & (cutName_cut_.find(*it)->second);
1635     if( c->filled )
1636     {
1637     c->histo1.Fill( c->value );
1638     if( passedAllPreviousCuts(c->variableName) ) c->histo2.Fill( c->value );
1639     if( passedAllOtherSameAndLowerLevelCuts(c->variableName) ) c->histo3.Fill( c->value );
1640     if( passedAllOtherCuts(c->variableName) ) c->histo4.Fill( c->value );
1641     if( passedCut("all") ) c->histo5.Fill( c->value );
1642     }
1643     }
1644     return ret;
1645     }
1646    
1647     bool baseClass::writeCutHistos()
1648     {
1649     bool ret = true;
1650     for (vector<string>::iterator it = orderedCutNames_.begin();
1651     it != orderedCutNames_.end(); it++)
1652     {
1653     cut * c = & (cutName_cut_.find(*it)->second);
1654     c->histo1.Write();
1655     c->histo2.Write();
1656     c->histo3.Write();
1657     c->histo4.Write();
1658     c->histo5.Write();
1659     }
1660    
1661     // Any failure mode to implement?
1662     return ret;
1663     }
1664    
1665     bool baseClass::updateCutEffic()
1666     {
1667     bool ret = true;
1668     for (vector<string>::iterator it = orderedCutNames_.begin();
1669     it != orderedCutNames_.end(); it++)
1670     {
1671     cut * c = & (cutName_cut_.find(*it)->second);
1672     if( passedAllPreviousCuts(c->variableName) )
1673     {
1674     c->nEvtInput++;
1675     if( passedCut(c->variableName) ) c->nEvtPassed++;
1676     }
1677     }
1678     return ret;
1679     }
1680    
1681    
1682     bool baseClass::writeCutEfficFile()
1683     {
1684    
1685     bool ret = true;
1686    
1687     // Set bin labels for event counter histogram
1688     int bincounter=1;
1689     eventcuts_->GetXaxis()->SetBinLabel(bincounter,"NoCuts");
1690     ++bincounter;
1691     if (skimWasMade_)
1692     {
1693     eventcuts_->GetXaxis()->SetBinLabel(bincounter,"Skim");
1694     ++bincounter;
1695     }
1696     for (int i=0;i<orderedCutNames_.size();++i)
1697     {
1698     eventcuts_->GetXaxis()->SetBinLabel(bincounter,orderedCutNames_[i].c_str());
1699     ++bincounter;
1700     }
1701    
1702     bincounter=1;
1703    
1704     int nEntRoottuple = fChain->GetEntriesFast();
1705     int nEntTot = (skimWasMade_ ? NBeforeSkim_ : nEntRoottuple );
1706     string cutEfficFile = *cutEfficFile_ + ".dat";
1707     ofstream os(cutEfficFile.c_str());
1708    
1709     os << "################################## Preliminary Cut Values ###################################################################\n"
1710     << "########################### variableName value1 value2 value3 value4 level\n"
1711     << preCutInfo_.str();
1712    
1713     int cutIdPed=0;
1714     os.precision(4);
1715     os << "################################## Cuts #####################################################################################\n"
1716     <<"#id variableName min1 max1 min2 max2 level N Npass EffRel errEffRel EffAbs errEffAbs"<<endl
1717     << fixed
1718     << setw(3) << cutIdPed
1719     << setw(25) << "nocut"
1720     << setprecision(4)
1721     << setw(15) << "-"
1722     << setw(15) << "-"
1723     << setw(15) << "-"
1724     << setw(15) << "-"
1725     << setw(15) << "-"
1726     << setw(15) << nEntTot
1727     << setw(15) << nEntTot
1728     << setprecision(11)
1729     << setw(15) << "1.00000000000"
1730     << setw(15) << "0.00000000000"
1731     << setw(15) << "1.00000000000"
1732     << setw(15) << "0.00000000000"
1733     << endl;
1734    
1735     double effRel;
1736     double effRelErr;
1737     double effAbs;
1738     double effAbsErr;
1739    
1740     eventcuts_->SetBinContent(bincounter,nEntTot);
1741     if (optimizeName_cut_.size())
1742     h_optimizer_->SetBinContent(0, nEntTot);
1743    
1744     if(skimWasMade_)
1745     {
1746     ++bincounter;
1747     eventcuts_->SetBinContent(bincounter, nEntRoottuple);
1748     effRel = (double) nEntRoottuple / (double) NBeforeSkim_;
1749     effRelErr = sqrt( (double) effRel * (1.0 - (double) effRel) / (double) NBeforeSkim_ );
1750     effAbs = effRel;
1751     effAbsErr = effRelErr;
1752     os << fixed
1753     << setw(3) << ++cutIdPed
1754     << setw(25) << "skim"
1755     << setprecision(4)
1756     << setw(15) << "-"
1757     << setw(15) << "-"
1758     << setw(15) << "-"
1759     << setw(15) << "-"
1760     << setw(15) << "-"
1761     << setw(15) << NBeforeSkim_
1762     << setw(15) << nEntRoottuple
1763     << setprecision(11)
1764     << setw(15) << effRel
1765     << setw(15) << effRelErr
1766     << setw(15) << effAbs
1767     << setw(15) << effAbsErr
1768     << endl;
1769     }
1770     for (vector<string>::iterator it = orderedCutNames_.begin();
1771     it != orderedCutNames_.end(); it++)
1772     {
1773     cut * c = & (cutName_cut_.find(*it)->second);
1774     ++bincounter;
1775     eventcuts_->SetBinContent(bincounter, c->nEvtPassed);
1776     effRel = (double) c->nEvtPassed / (double) c->nEvtInput;
1777     effRelErr = sqrt( (double) effRel * (1.0 - (double) effRel) / (double) c->nEvtInput );
1778     effAbs = (double) c->nEvtPassed / (double) nEntTot;
1779     effAbsErr = sqrt( (double) effAbs * (1.0 - (double) effAbs) / (double) nEntTot );
1780    
1781     std::stringstream ssm1, ssM1, ssm2,ssM2;
1782     ssm1 << fixed << setprecision(4) << c->minValue1;
1783     ssM1 << fixed << setprecision(4) << c->maxValue1;
1784     if(c->minValue2 == -9999999)
1785     {
1786     ssm2 << "-inf";
1787     }
1788     else
1789     {
1790     ssm2 << fixed << setprecision(4) << c->minValue2;
1791     }
1792     if(c->maxValue2 == 9999999)
1793     {
1794     ssM2 << "+inf";
1795     }
1796     else
1797     {
1798     ssM2 << fixed << setprecision(4) << c->maxValue2;
1799     }
1800     os << setw(3) << cutIdPed+c->id
1801     << setw(25) << c->variableName
1802     << setprecision(4)
1803     << fixed
1804     << setw(15) << ( ( c->minValue1 == -9999999.0 ) ? "-inf" : ssm1.str() )
1805     << setw(15) << ( ( c->maxValue1 == 9999999.0 ) ? "+inf" : ssM1.str() )
1806     << setw(15) << ( ( c->minValue2 > c->maxValue2 ) ? "-" : ssm2.str() )
1807     << setw(15) << ( ( c->minValue2 > c->maxValue2 ) ? "-" : ssM2.str() )
1808     << setw(15) << c->level_int
1809     << setw(15) << c->nEvtInput
1810     << setw(15) << c->nEvtPassed
1811     << setprecision(11)
1812     << setw(15) << effRel
1813     << setw(15) << effRelErr
1814     << setw(15) << effAbs
1815     << setw(15) << effAbsErr
1816     << endl;
1817     }
1818    
1819     // Write optimization histograms
1820     if (optimizeName_cut_.size())
1821     {
1822     gDirectory->mkdir("Optimizer");
1823     gDirectory->cd("Optimizer");
1824     h_optimizer_->Write();
1825     for (int i=0;i<optimizeName_cut_.size();++i)
1826     {
1827     stringstream x;
1828     x<<"Cut"<<i<<"_"<<optimizeName_cut_[i].variableName;
1829     if (optimizeName_cut_[i].testgreater==true)
1830     x<<"_gt_";
1831     else
1832     x<<"_lt_";
1833     x<<optimizeName_cut_[i].minvalue<<"_to_"<<optimizeName_cut_[i].maxvalue;
1834     TObjString test(x.str().c_str());
1835     test.Write();
1836     }
1837     gDirectory->cd("..");
1838     }
1839    
1840     eventcuts_->Write(); // write here, since WriteCutHistos is called before WriteCutEfficFile
1841     // Any failure mode to implement?
1842     return ret;
1843     } // writeCutEffFile
1844    
1845    
1846     bool baseClass::sortCuts(const cut& X, const cut& Y)
1847     {
1848     return X.id < Y.id;
1849     }
1850    
1851    
1852     vector<string> baseClass::split(const string& s)
1853     {
1854     vector<string> ret;
1855     string::size_type i =0;
1856     while (i != s.size()){
1857     while (i != s.size() && isspace(s[i]))
1858     ++i;
1859     string::size_type j = i;
1860     while (j != s.size() && !isspace(s[j]))
1861     ++j;
1862     if (i != j){
1863     ret.push_back(s.substr(i, j -i));
1864     i = j;
1865     }
1866     }
1867     return ret;
1868     }
1869    
1870     double baseClass::decodeCutValue(const string& s)
1871     {
1872     //STDOUT("s = "<<s);
1873     double ret;
1874     if( s == "inf" || s == "+inf" )
1875     {
1876     ret = 9999999;
1877     }
1878     else if ( s == "-inf" || s == "-" )
1879     {
1880     ret = -9999999;
1881     }
1882     else
1883     {
1884     ret = atof( s.c_str() );
1885     }
1886     return ret;
1887     }
1888    
1889     int baseClass::getGlobalInfoNstart(char *pName)
1890     {
1891     int NBeforeSkim = 0;
1892     STDOUT(pName<<" "<< NBeforeSkim)
1893     TFile *f = new TFile(pName);
1894     string s1 = "LJFilter/EventCount/EventCounter";
1895     string s2 = "LJFilterPAT/EventCount/EventCounter";
1896     TH1I* hCount1 = (TH1I*)f->Get(s1.c_str());
1897     TH1I* hCount2 = (TH1I*)f->Get(s2.c_str());
1898     if( !hCount1 && !hCount2 )
1899     {
1900     STDOUT("Skim filter histogram(s) not found. Will assume skim was not made for ALL files.");
1901     skimWasMade_ = false;
1902     return NBeforeSkim;
1903     }
1904    
1905     if (hCount1) NBeforeSkim = (int)hCount1->GetBinContent(1);
1906     else NBeforeSkim = (int)hCount2->GetBinContent(1);
1907    
1908     // STDOUT(pName<<" "<< NBeforeSkim)
1909     f->Close();
1910    
1911     return NBeforeSkim;
1912     }
1913    
1914