ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/ElectronEnergyRegression.cc
Revision: 1.1
Committed: Tue Oct 9 13:57:28 2012 UTC (12 years, 6 months ago) by sixie
Content type: text/plain
Branch: MAIN
Log Message:
add electron energy regression evaluator

File Contents

# User Rev Content
1 sixie 1.1 #include "MitPhysics/Utils/interface/ElectronEnergyRegression.h"
2     #include <cmath>
3     #include <cassert>
4     #include <iostream>
5    
6     ClassImp(mithep::ElectronEnergyRegression)
7    
8     using namespace mithep;
9    
10    
11     ElectronEnergyRegression::ElectronEnergyRegression() :
12     fIsInitialized(kFALSE),
13     fVersionType(kNoTrkVar),
14     forestCorrection_eb(0),
15     forestCorrection_ee(0),
16     forestUncertainty_eb(0),
17     forestUncertainty_ee(0) {
18     }
19    
20     ElectronEnergyRegression::~ElectronEnergyRegression() {}
21     // Destructor does nothing
22    
23    
24     void ElectronEnergyRegression::initialize(std::string weightsFile,
25     ElectronEnergyRegression::ElectronEnergyRegressionType type) {
26    
27     // Loading forest object according to different versions
28     TFile file(weightsFile.c_str());
29    
30     if (type == kNoTrkVar || type == kWithTrkVar) {
31     forestCorrection_eb = (GBRForest*) file.Get("EBCorrection");
32     forestCorrection_ee = (GBRForest*) file.Get("EECorrection");
33     forestUncertainty_eb = (GBRForest*) file.Get("EBUncertainty");
34     forestUncertainty_ee = (GBRForest*) file.Get("EEUncertainty");
35    
36     // Just checking
37     assert(forestCorrection_eb);
38     assert(forestCorrection_ee);
39     assert(forestUncertainty_eb);
40     assert(forestUncertainty_ee);
41     }
42    
43     // Updating type and marking as initialized
44     fVersionType = type;
45     fIsInitialized = kTRUE;
46     }
47    
48    
49     double ElectronEnergyRegression::regressionValueNoTrkVar(
50     double SCRawEnergy,
51     double scEta,
52     double scPhi,
53     double R9,
54     double etawidth,
55     double phiwidth,
56     double NClusters,
57     double HoE,
58     double rho,
59     double vertices,
60     double EtaSeed,
61     double PhiSeed,
62     double ESeed,
63     double E3x3Seed,
64     double E5x5Seed,
65     double see,
66     double spp,
67     double sep,
68     double EMaxSeed,
69     double E2ndSeed,
70     double ETopSeed,
71     double EBottomSeed,
72     double ELeftSeed,
73     double ERightSeed,
74     double E2x5MaxSeed,
75     double E2x5TopSeed,
76     double E2x5BottomSeed,
77     double E2x5LeftSeed,
78     double E2x5RightSeed,
79     double IEtaSeed,
80     double IPhiSeed,
81     double EtaCrySeed,
82     double PhiCrySeed,
83     double PreShowerOverRaw,
84     bool printDebug)
85     {
86     // Checking if instance has been initialized
87     if (fIsInitialized == kFALSE) {
88     printf("ElectronEnergyRegression instance not initialized !!!");
89     return 0;
90     }
91    
92     // Checking if type is correct
93     if (!(fVersionType == kNoTrkVar)) {
94     std::cout << "Error: Regression VersionType " << fVersionType << " is not supported to use function regressionValueNoTrkVar.\n";
95     return 0;
96     }
97    
98     // Now applying regression according to version and (endcap/barrel)
99     float *vals = (fabs(scEta) <= 1.479) ? new float[38] : new float[31];
100     if (fabs(scEta) <= 1.479) { // Barrel
101     vals[0] = SCRawEnergy;
102     vals[1] = scEta;
103     vals[2] = scPhi;
104     vals[3] = R9;
105     vals[4] = E5x5Seed/SCRawEnergy;
106     vals[5] = etawidth;
107     vals[6] = phiwidth;
108     vals[7] = NClusters;
109     vals[8] = HoE;
110     vals[9] = rho;
111     vals[10] = vertices;
112     vals[11] = EtaSeed - scEta;
113     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
114     vals[13] = ESeed/SCRawEnergy;
115     vals[14] = E3x3Seed/ESeed;
116     vals[15] = E5x5Seed/ESeed;
117     vals[16] = see;
118     vals[17] = spp;
119     vals[18] = sep;
120     vals[19] = EMaxSeed/ESeed;
121     vals[20] = E2ndSeed/ESeed;
122     vals[21] = ETopSeed/ESeed;
123     vals[22] = EBottomSeed/ESeed;
124     vals[23] = ELeftSeed/ESeed;
125     vals[24] = ERightSeed/ESeed;
126     vals[25] = E2x5MaxSeed/ESeed;
127     vals[26] = E2x5TopSeed/ESeed;
128     vals[27] = E2x5BottomSeed/ESeed;
129     vals[28] = E2x5LeftSeed/ESeed;
130     vals[29] = E2x5RightSeed/ESeed;
131     vals[30] = IEtaSeed;
132     vals[31] = IPhiSeed;
133     vals[32] = ((int) IEtaSeed)%5;
134     vals[33] = ((int) IPhiSeed)%2;
135     vals[34] = (abs(IEtaSeed)<=25)*(((int)IEtaSeed)%25) + (abs(IEtaSeed)>25)*(((int) (IEtaSeed-25*abs(IEtaSeed)/IEtaSeed))%20);
136     vals[35] = ((int) IPhiSeed)%20;
137     vals[36] = EtaCrySeed;
138     vals[37] = PhiCrySeed;
139     }
140     else { // Endcap
141     vals[0] = SCRawEnergy;
142     vals[1] = scEta;
143     vals[2] = scPhi;
144     vals[3] = R9;
145     vals[4] = E5x5Seed/SCRawEnergy;
146     vals[5] = etawidth;
147     vals[6] = phiwidth;
148     vals[7] = NClusters;
149     vals[8] = HoE;
150     vals[9] = rho;
151     vals[10] = vertices;
152     vals[11] = EtaSeed - scEta;
153     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
154     vals[13] = ESeed/SCRawEnergy;
155     vals[14] = E3x3Seed/ESeed;
156     vals[15] = E5x5Seed/ESeed;
157     vals[16] = see;
158     vals[17] = spp;
159     vals[18] = sep;
160     vals[19] = EMaxSeed/ESeed;
161     vals[20] = E2ndSeed/ESeed;
162     vals[21] = ETopSeed/ESeed;
163     vals[22] = EBottomSeed/ESeed;
164     vals[23] = ELeftSeed/ESeed;
165     vals[24] = ERightSeed/ESeed;
166     vals[25] = E2x5MaxSeed/ESeed;
167     vals[26] = E2x5TopSeed/ESeed;
168     vals[27] = E2x5BottomSeed/ESeed;
169     vals[28] = E2x5LeftSeed/ESeed;
170     vals[29] = E2x5RightSeed/ESeed;
171     vals[30] = PreShowerOverRaw;
172     }
173    
174     // Now evaluating the regression
175     double regressionResult = 0;
176     Int_t BinIndex = -1;
177    
178     if (fVersionType == kNoTrkVar) {
179     if (fabs(scEta) <= 1.479) {
180     regressionResult = SCRawEnergy * forestCorrection_eb->GetResponse(vals);
181     BinIndex = 0;
182     }
183     else {
184     regressionResult = (SCRawEnergy*(1+PreShowerOverRaw)) * forestCorrection_ee->GetResponse(vals);
185     BinIndex = 1;
186     }
187     }
188    
189     //print debug
190     if (printDebug) {
191     if ( fabs(scEta) <= 1.479) {
192     std::cout << "Barrel :";
193     for (uint v=0; v < 38; ++v) std::cout << vals[v] << ", ";
194     std::cout << "\n";
195     }
196     else {
197     std::cout << "Endcap :";
198     for (uint v=0; v < 32; ++v) std::cout << vals[v] << ", ";
199     std::cout << "\n";
200     }
201     std::cout << "BinIndex : " << BinIndex << "\n";
202     std::cout << "SCRawEnergy = " << SCRawEnergy << " : PreShowerOverRaw = " << PreShowerOverRaw << std::endl;
203     std::cout << "regression energy = " << regressionResult << std::endl;
204     }
205    
206    
207     // Cleaning up and returning
208     delete[] vals;
209     return regressionResult;
210     }
211    
212     double ElectronEnergyRegression::regressionUncertaintyNoTrkVar(
213     double SCRawEnergy,
214     double scEta,
215     double scPhi,
216     double R9,
217     double etawidth,
218     double phiwidth,
219     double NClusters,
220     double HoE,
221     double rho,
222     double vertices,
223     double EtaSeed,
224     double PhiSeed,
225     double ESeed,
226     double E3x3Seed,
227     double E5x5Seed,
228     double see,
229     double spp,
230     double sep,
231     double EMaxSeed,
232     double E2ndSeed,
233     double ETopSeed,
234     double EBottomSeed,
235     double ELeftSeed,
236     double ERightSeed,
237     double E2x5MaxSeed,
238     double E2x5TopSeed,
239     double E2x5BottomSeed,
240     double E2x5LeftSeed,
241     double E2x5RightSeed,
242     double IEtaSeed,
243     double IPhiSeed,
244     double EtaCrySeed,
245     double PhiCrySeed,
246     double PreShowerOverRaw,
247     bool printDebug)
248     {
249     // Checking if instance has been initialized
250     if (fIsInitialized == kFALSE) {
251     printf("ElectronEnergyRegression instance not initialized !!!");
252     return 0;
253     }
254    
255     // Checking if type is correct
256     if (!(fVersionType == kNoTrkVar)) {
257     std::cout << "Error: Regression VersionType " << fVersionType << " is not supported to use function regressionValueNoTrkVar.\n";
258     return 0;
259     }
260    
261     // Now applying regression according to version and (endcap/barrel)
262     float *vals = (fabs(scEta) <= 1.479) ? new float[38] : new float[31];
263     if (fabs(scEta) <= 1.479) { // Barrel
264     vals[0] = SCRawEnergy;
265     vals[1] = scEta;
266     vals[2] = scPhi;
267     vals[3] = R9;
268     vals[4] = E5x5Seed/SCRawEnergy;
269     vals[5] = etawidth;
270     vals[6] = phiwidth;
271     vals[7] = NClusters;
272     vals[8] = HoE;
273     vals[9] = rho;
274     vals[10] = vertices;
275     vals[11] = EtaSeed - scEta;
276     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
277     vals[13] = ESeed/SCRawEnergy;
278     vals[14] = E3x3Seed/ESeed;
279     vals[15] = E5x5Seed/ESeed;
280     vals[16] = see;
281     vals[17] = spp;
282     vals[18] = sep;
283     vals[19] = EMaxSeed/ESeed;
284     vals[20] = E2ndSeed/ESeed;
285     vals[21] = ETopSeed/ESeed;
286     vals[22] = EBottomSeed/ESeed;
287     vals[23] = ELeftSeed/ESeed;
288     vals[24] = ERightSeed/ESeed;
289     vals[25] = E2x5MaxSeed/ESeed;
290     vals[26] = E2x5TopSeed/ESeed;
291     vals[27] = E2x5BottomSeed/ESeed;
292     vals[28] = E2x5LeftSeed/ESeed;
293     vals[29] = E2x5RightSeed/ESeed;
294     vals[30] = IEtaSeed;
295     vals[31] = IPhiSeed;
296     vals[32] = ((int) IEtaSeed)%5;
297     vals[33] = ((int) IPhiSeed)%2;
298     vals[34] = (abs(IEtaSeed)<=25)*(((int)IEtaSeed)%25) + (abs(IEtaSeed)>25)*(((int) (IEtaSeed-25*abs(IEtaSeed)/IEtaSeed))%20);
299     vals[35] = ((int) IPhiSeed)%20;
300     vals[36] = EtaCrySeed;
301     vals[37] = PhiCrySeed;
302     }
303     else { // Endcap
304     vals[0] = SCRawEnergy;
305     vals[1] = scEta;
306     vals[2] = scPhi;
307     vals[3] = R9;
308     vals[4] = E5x5Seed/SCRawEnergy;
309     vals[5] = etawidth;
310     vals[6] = phiwidth;
311     vals[7] = NClusters;
312     vals[8] = HoE;
313     vals[9] = rho;
314     vals[10] = vertices;
315     vals[11] = EtaSeed - scEta;
316     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
317     vals[13] = ESeed/SCRawEnergy;
318     vals[14] = E3x3Seed/ESeed;
319     vals[15] = E5x5Seed/ESeed;
320     vals[16] = see;
321     vals[17] = spp;
322     vals[18] = sep;
323     vals[19] = EMaxSeed/ESeed;
324     vals[20] = E2ndSeed/ESeed;
325     vals[21] = ETopSeed/ESeed;
326     vals[22] = EBottomSeed/ESeed;
327     vals[23] = ELeftSeed/ESeed;
328     vals[24] = ERightSeed/ESeed;
329     vals[25] = E2x5MaxSeed/ESeed;
330     vals[26] = E2x5TopSeed/ESeed;
331     vals[27] = E2x5BottomSeed/ESeed;
332     vals[28] = E2x5LeftSeed/ESeed;
333     vals[29] = E2x5RightSeed/ESeed;
334     vals[30] = PreShowerOverRaw;
335     }
336    
337     // Now evaluating the regression
338     double regressionResult = 0;
339     Int_t BinIndex = -1;
340    
341     if (fVersionType == kNoTrkVar) {
342     if (fabs(scEta) <= 1.479) {
343     regressionResult = SCRawEnergy * forestUncertainty_eb->GetResponse(vals);
344     BinIndex = 0;
345     }
346     else {
347     regressionResult = (SCRawEnergy*(1+PreShowerOverRaw)) * forestUncertainty_ee->GetResponse(vals);
348     BinIndex = 1;
349     }
350     }
351    
352     //print debug
353     if (printDebug) {
354     if (fabs(scEta) <= 1.479) {
355     std::cout << "Barrel :";
356     for (uint v=0; v < 38; ++v) std::cout << vals[v] << ", ";
357     std::cout << "\n";
358     }
359     else {
360     std::cout << "Endcap :";
361     for (uint v=0; v < 32; ++v) std::cout << vals[v] << ", ";
362     std::cout << "\n";
363     }
364     std::cout << "BinIndex : " << BinIndex << "\n";
365     std::cout << "SCRawEnergy = " << SCRawEnergy << " : PreShowerOverRaw = " << PreShowerOverRaw << std::endl;
366     std::cout << "regression energy uncertainty = " << regressionResult << std::endl;
367     }
368    
369    
370     // Cleaning up and returning
371     delete[] vals;
372     return regressionResult;
373     }
374    
375     double ElectronEnergyRegression::regressionValueWithTrkVarV1(
376     double SCRawEnergy,
377     double scEta,
378     double scPhi,
379     double R9,
380     double etawidth,
381     double phiwidth,
382     double NClusters,
383     double HoE,
384     double rho,
385     double vertices,
386     double EtaSeed,
387     double PhiSeed,
388     double ESeed,
389     double E3x3Seed,
390     double E5x5Seed,
391     double see,
392     double spp,
393     double sep,
394     double EMaxSeed,
395     double E2ndSeed,
396     double ETopSeed,
397     double EBottomSeed,
398     double ELeftSeed,
399     double ERightSeed,
400     double E2x5MaxSeed,
401     double E2x5TopSeed,
402     double E2x5BottomSeed,
403     double E2x5LeftSeed,
404     double E2x5RightSeed,
405     double IEtaSeed,
406     double IPhiSeed,
407     double EtaCrySeed,
408     double PhiCrySeed,
409     double PreShowerOverRaw,
410     double GsfTrackPIn,
411     double fbrem,
412     double Charge,
413     double EoP,
414     double TrackMomentumError,
415     bool printDebug)
416     {
417     // Checking if instance has been initialized
418     if (fIsInitialized == kFALSE) {
419     printf("ElectronEnergyRegression instance not initialized !!!");
420     return 0;
421     }
422    
423     // Checking if fVersionType is correct
424     assert(fVersionType == kWithTrkVar);
425    
426     float *vals = (fabs(scEta) <= 1.479) ? new float[43] : new float[36];
427     if (fabs(scEta) <= 1.479) { // Barrel
428     vals[0] = SCRawEnergy;
429     vals[1] = scEta;
430     vals[2] = scPhi;
431     vals[3] = R9;
432     vals[4] = E5x5Seed/SCRawEnergy;
433     vals[5] = etawidth;
434     vals[6] = phiwidth;
435     vals[7] = NClusters;
436     vals[8] = HoE;
437     vals[9] = rho;
438     vals[10] = vertices;
439     vals[11] = EtaSeed - scEta;
440     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
441     vals[13] = ESeed/SCRawEnergy;
442     vals[14] = E3x3Seed/ESeed;
443     vals[15] = E5x5Seed/ESeed;
444     vals[16] = see;
445     vals[17] = spp;
446     vals[18] = sep;
447     vals[19] = EMaxSeed/ESeed;
448     vals[20] = E2ndSeed/ESeed;
449     vals[21] = ETopSeed/ESeed;
450     vals[22] = EBottomSeed/ESeed;
451     vals[23] = ELeftSeed/ESeed;
452     vals[24] = ERightSeed/ESeed;
453     vals[25] = E2x5MaxSeed/ESeed;
454     vals[26] = E2x5TopSeed/ESeed;
455     vals[27] = E2x5BottomSeed/ESeed;
456     vals[28] = E2x5LeftSeed/ESeed;
457     vals[29] = E2x5RightSeed/ESeed;
458     vals[30] = GsfTrackPIn;
459     vals[31] = fbrem;
460     vals[32] = Charge;
461     vals[33] = EoP;
462     vals[34] = TrackMomentumError/GsfTrackPIn;
463     vals[35] = IEtaSeed;
464     vals[36] = IPhiSeed;
465     vals[37] = ((int) IEtaSeed)%5;
466     vals[38] = ((int) IPhiSeed)%2;
467     vals[39] = (abs(IEtaSeed)<=25)*(((int)IEtaSeed)%25) + (abs(IEtaSeed)>25)*(((int) (IEtaSeed-25*abs(IEtaSeed)/IEtaSeed))%20);
468     vals[40] = ((int) IPhiSeed)%20;
469     vals[41] = EtaCrySeed;
470     vals[42] = PhiCrySeed;
471     }
472    
473     else { // Endcap
474     vals[0] = SCRawEnergy;
475     vals[1] = scEta;
476     vals[2] = scPhi;
477     vals[3] = R9;
478     vals[4] = E5x5Seed/SCRawEnergy;
479     vals[5] = etawidth;
480     vals[6] = phiwidth;
481     vals[7] = NClusters;
482     vals[8] = HoE;
483     vals[9] = rho;
484     vals[10] = vertices;
485     vals[11] = EtaSeed - scEta;
486     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
487     vals[13] = ESeed/SCRawEnergy;
488     vals[14] = E3x3Seed/ESeed;
489     vals[15] = E5x5Seed/ESeed;
490     vals[16] = see;
491     vals[17] = spp;
492     vals[18] = sep;
493     vals[19] = EMaxSeed/ESeed;
494     vals[20] = E2ndSeed/ESeed;
495     vals[21] = ETopSeed/ESeed;
496     vals[22] = EBottomSeed/ESeed;
497     vals[23] = ELeftSeed/ESeed;
498     vals[24] = ERightSeed/ESeed;
499     vals[25] = E2x5MaxSeed/ESeed;
500     vals[26] = E2x5TopSeed/ESeed;
501     vals[27] = E2x5BottomSeed/ESeed;
502     vals[28] = E2x5LeftSeed/ESeed;
503     vals[29] = E2x5RightSeed/ESeed;
504     vals[30] = GsfTrackPIn;
505     vals[31] = fbrem;
506     vals[32] = Charge;
507     vals[33] = EoP;
508     vals[34] = TrackMomentumError/GsfTrackPIn;
509     vals[35] = PreShowerOverRaw;
510     }
511    
512     // Now evaluating the regression
513     double regressionResult = 0;
514    
515     if (fVersionType == kWithTrkVar) {
516     if (fabs(scEta) <= 1.479) regressionResult = SCRawEnergy * forestCorrection_eb->GetResponse(vals);
517     else regressionResult = (SCRawEnergy*(1+PreShowerOverRaw)) * forestCorrection_ee->GetResponse(vals);
518     }
519    
520    
521     //print debug
522     if (printDebug) {
523     if (scEta <= 1.479) {
524     std::cout << "Barrel :";
525     for (uint v=0; v < 43; ++v) std::cout << vals[v] << ", ";
526     std::cout << "\n";
527     }
528     else {
529     std::cout << "Endcap :";
530     for (uint v=0; v < 36; ++v) std::cout << vals[v] << ", ";
531     std::cout << "\n";
532     }
533     std::cout << "SCRawEnergy = " << SCRawEnergy << " : PreShowerOverRaw = " << PreShowerOverRaw << std::endl;
534     std::cout << "regression energy = " << regressionResult << std::endl;
535     }
536    
537     // Cleaning up and returning
538     delete[] vals;
539     return regressionResult;
540     }
541    
542    
543    
544    
545     double ElectronEnergyRegression::regressionUncertaintyWithTrkVarV1(
546     double SCRawEnergy,
547     double scEta,
548     double scPhi,
549     double R9,
550     double etawidth,
551     double phiwidth,
552     double NClusters,
553     double HoE,
554     double rho,
555     double vertices,
556     double EtaSeed,
557     double PhiSeed,
558     double ESeed,
559     double E3x3Seed,
560     double E5x5Seed,
561     double see,
562     double spp,
563     double sep,
564     double EMaxSeed,
565     double E2ndSeed,
566     double ETopSeed,
567     double EBottomSeed,
568     double ELeftSeed,
569     double ERightSeed,
570     double E2x5MaxSeed,
571     double E2x5TopSeed,
572     double E2x5BottomSeed,
573     double E2x5LeftSeed,
574     double E2x5RightSeed,
575     double IEtaSeed,
576     double IPhiSeed,
577     double EtaCrySeed,
578     double PhiCrySeed,
579     double PreShowerOverRaw,
580     double GsfTrackPIn,
581     double fbrem,
582     double Charge,
583     double EoP,
584     double TrackMomentumError,
585     bool printDebug)
586     {
587     // Checking if instance has been initialized
588     if (fIsInitialized == kFALSE) {
589     printf("ElectronEnergyRegression instance not initialized !!!");
590     return 0;
591     }
592    
593     // Checking if fVersionType is correct
594     assert(fVersionType == kWithTrkVar);
595    
596     float *vals = (fabs(scEta) <= 1.479) ? new float[43] : new float[36];
597     if (fabs(scEta) <= 1.479) { // Barrel
598     vals[0] = SCRawEnergy;
599     vals[1] = scEta;
600     vals[2] = scPhi;
601     vals[3] = R9;
602     vals[4] = E5x5Seed/SCRawEnergy;
603     vals[5] = etawidth;
604     vals[6] = phiwidth;
605     vals[7] = NClusters;
606     vals[8] = HoE;
607     vals[9] = rho;
608     vals[10] = vertices;
609     vals[11] = EtaSeed - scEta;
610     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
611     vals[13] = ESeed/SCRawEnergy;
612     vals[14] = E3x3Seed/ESeed;
613     vals[15] = E5x5Seed/ESeed;
614     vals[16] = see;
615     vals[17] = spp;
616     vals[18] = sep;
617     vals[19] = EMaxSeed/ESeed;
618     vals[20] = E2ndSeed/ESeed;
619     vals[21] = ETopSeed/ESeed;
620     vals[22] = EBottomSeed/ESeed;
621     vals[23] = ELeftSeed/ESeed;
622     vals[24] = ERightSeed/ESeed;
623     vals[25] = E2x5MaxSeed/ESeed;
624     vals[26] = E2x5TopSeed/ESeed;
625     vals[27] = E2x5BottomSeed/ESeed;
626     vals[28] = E2x5LeftSeed/ESeed;
627     vals[29] = E2x5RightSeed/ESeed;
628     vals[30] = GsfTrackPIn;
629     vals[31] = fbrem;
630     vals[32] = Charge;
631     vals[33] = EoP;
632     vals[34] = TrackMomentumError/GsfTrackPIn;
633     vals[35] = IEtaSeed;
634     vals[36] = IPhiSeed;
635     vals[37] = ((int) IEtaSeed)%5;
636     vals[38] = ((int) IPhiSeed)%2;
637     vals[39] = (abs(IEtaSeed)<=25)*(((int)IEtaSeed)%25) + (abs(IEtaSeed)>25)*(((int) (IEtaSeed-25*abs(IEtaSeed)/IEtaSeed))%20);
638     vals[40] = ((int) IPhiSeed)%20;
639     vals[41] = EtaCrySeed;
640     vals[42] = PhiCrySeed;
641     }
642    
643     else { // Endcap
644     vals[0] = SCRawEnergy;
645     vals[1] = scEta;
646     vals[2] = scPhi;
647     vals[3] = R9;
648     vals[4] = E5x5Seed/SCRawEnergy;
649     vals[5] = etawidth;
650     vals[6] = phiwidth;
651     vals[7] = NClusters;
652     vals[8] = HoE;
653     vals[9] = rho;
654     vals[10] = vertices;
655     vals[11] = EtaSeed - scEta;
656     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
657     vals[13] = ESeed/SCRawEnergy;
658     vals[14] = E3x3Seed/ESeed;
659     vals[15] = E5x5Seed/ESeed;
660     vals[16] = see;
661     vals[17] = spp;
662     vals[18] = sep;
663     vals[19] = EMaxSeed/ESeed;
664     vals[20] = E2ndSeed/ESeed;
665     vals[21] = ETopSeed/ESeed;
666     vals[22] = EBottomSeed/ESeed;
667     vals[23] = ELeftSeed/ESeed;
668     vals[24] = ERightSeed/ESeed;
669     vals[25] = E2x5MaxSeed/ESeed;
670     vals[26] = E2x5TopSeed/ESeed;
671     vals[27] = E2x5BottomSeed/ESeed;
672     vals[28] = E2x5LeftSeed/ESeed;
673     vals[29] = E2x5RightSeed/ESeed;
674     vals[30] = GsfTrackPIn;
675     vals[31] = fbrem;
676     vals[32] = Charge;
677     vals[33] = EoP;
678     vals[34] = TrackMomentumError/GsfTrackPIn;
679     vals[35] = PreShowerOverRaw;
680     }
681    
682     // Now evaluating the regression
683     double regressionResult = 0;
684    
685     if (fVersionType == kWithTrkVar) {
686     if (fabs(scEta) <= 1.479) regressionResult = SCRawEnergy * forestUncertainty_eb->GetResponse(vals);
687     else regressionResult = (SCRawEnergy*(1+PreShowerOverRaw)) * forestUncertainty_ee->GetResponse(vals);
688     }
689    
690     //print debug
691     if (printDebug) {
692     if (scEta <= 1.479) {
693     std::cout << "Barrel :";
694     for (uint v=0; v < 43; ++v) std::cout << vals[v] << ", ";
695     std::cout << "\n";
696     }
697     else {
698     std::cout << "Endcap :";
699     for (uint v=0; v < 36; ++v) std::cout << vals[v] << ", ";
700     std::cout << "\n";
701     }
702     std::cout << " SCRawEnergy = " << SCRawEnergy << " : PreShowerOverRaw = " << PreShowerOverRaw << std::endl;
703     std::cout << "regression energy uncertainty = " << regressionResult << std::endl;
704     }
705    
706    
707     // Cleaning up and returning
708     delete[] vals;
709     return regressionResult;
710     }
711    
712    
713    
714    
715     double ElectronEnergyRegression::regressionValueWithTrkVarV2( std::vector<double> &inputvars,
716     bool printDebug)
717     {
718     // Checking if instance has been initialized
719     if (fIsInitialized == kFALSE) {
720     printf("ElectronEnergyRegression instance not initialized !!!");
721     return 0;
722     }
723    
724     //Check that inputvars vector has the right number of inputs
725     assert(inputvars.size() == 46);
726    
727     // Checking if fVersionType is correct
728     assert(fVersionType == kWithTrkVar);
729    
730     //assign variables from inputvars array to named variables
731     double SCRawEnergy = inputvars[0];
732     double scEta = inputvars[1];
733     double scPhi = inputvars[2];
734     double R9 = inputvars[3];
735     double etawidth = inputvars[4];
736     double phiwidth = inputvars[5];
737     double NClusters = inputvars[6];
738     double HoE = inputvars[7];
739     double rho = inputvars[8];
740     double vertices = inputvars[9];
741     double EtaSeed = inputvars[10];
742     double PhiSeed = inputvars[11];
743     double ESeed = inputvars[12];
744     double E3x3Seed = inputvars[13];
745     double E5x5Seed = inputvars[14];
746     double see = inputvars[15];
747     double spp = inputvars[16];
748     double sep = inputvars[17];
749     double EMaxSeed = inputvars[18];
750     double E2ndSeed = inputvars[19];
751     double ETopSeed = inputvars[20];
752     double EBottomSeed = inputvars[21];
753     double ELeftSeed = inputvars[22];
754     double ERightSeed = inputvars[23];
755     double E2x5MaxSeed = inputvars[24];
756     double E2x5TopSeed = inputvars[25];
757     double E2x5BottomSeed = inputvars[26];
758     double E2x5LeftSeed = inputvars[27];
759     double E2x5RightSeed = inputvars[28];
760     double IEtaSeed = inputvars[29];
761     double IPhiSeed = inputvars[30];
762     double EtaCrySeed = inputvars[31];
763     double PhiCrySeed = inputvars[32];
764     double PreShowerOverRaw = inputvars[33];
765     double GsfTrackPIn = inputvars[34];
766     double fbrem = inputvars[35];
767     double Charge = inputvars[36];
768     double EoP = inputvars[37];
769     double TrackMomentumError = inputvars[38];
770     double detaIn = inputvars[39];
771     double dphiIn = inputvars[40];
772     double detaCalo = inputvars[41];
773     double dphiCalo = inputvars[42];
774     double GsfTrackChiSqr = inputvars[43];
775     double KFTrackNLayers = inputvars[44];
776     double ElectronEnergyOverPout = inputvars[45];
777    
778    
779     float *vals = (fabs(scEta) <= 1.479) ? new float[50] : new float[43];
780     if (fabs(scEta) <= 1.479) { // Barrel
781     vals[0] = SCRawEnergy;
782     vals[1] = scEta;
783     vals[2] = scPhi;
784     vals[3] = R9;
785     vals[4] = E5x5Seed/SCRawEnergy;
786     vals[5] = etawidth;
787     vals[6] = phiwidth;
788     vals[7] = NClusters;
789     vals[8] = HoE;
790     vals[9] = rho;
791     vals[10] = vertices;
792     vals[11] = EtaSeed - scEta;
793     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
794     vals[13] = ESeed/SCRawEnergy;
795     vals[14] = E3x3Seed/ESeed;
796     vals[15] = E5x5Seed/ESeed;
797     vals[16] = see;
798     vals[17] = spp;
799     vals[18] = sep;
800     vals[19] = EMaxSeed/ESeed;
801     vals[20] = E2ndSeed/ESeed;
802     vals[21] = ETopSeed/ESeed;
803     vals[22] = EBottomSeed/ESeed;
804     vals[23] = ELeftSeed/ESeed;
805     vals[24] = ERightSeed/ESeed;
806     vals[25] = E2x5MaxSeed/ESeed;
807     vals[26] = E2x5TopSeed/ESeed;
808     vals[27] = E2x5BottomSeed/ESeed;
809     vals[28] = E2x5LeftSeed/ESeed;
810     vals[29] = E2x5RightSeed/ESeed;
811     vals[30] = GsfTrackPIn;
812     vals[31] = fbrem;
813     vals[32] = Charge;
814     vals[33] = EoP;
815     vals[34] = TrackMomentumError/GsfTrackPIn;
816     vals[35] = detaIn;
817     vals[36] = dphiIn;
818     vals[37] = detaCalo;
819     vals[38] = dphiCalo;
820     vals[39] = GsfTrackChiSqr;
821     vals[40] = KFTrackNLayers;
822     vals[41] = ElectronEnergyOverPout;
823     vals[42] = IEtaSeed;
824     vals[43] = IPhiSeed;
825     vals[44] = ((int) IEtaSeed)%5;
826     vals[45] = ((int) IPhiSeed)%2;
827     vals[46] = (abs(IEtaSeed)<=25)*(((int)IEtaSeed)%25) + (abs(IEtaSeed)>25)*(((int) (IEtaSeed-25*abs(IEtaSeed)/IEtaSeed))%20);
828     vals[47] = ((int) IPhiSeed)%20;
829     vals[48] = EtaCrySeed;
830     vals[49] = PhiCrySeed;
831     }
832    
833     else { // Endcap
834     vals[0] = SCRawEnergy;
835     vals[1] = scEta;
836     vals[2] = scPhi;
837     vals[3] = R9;
838     vals[4] = E5x5Seed/SCRawEnergy;
839     vals[5] = etawidth;
840     vals[6] = phiwidth;
841     vals[7] = NClusters;
842     vals[8] = HoE;
843     vals[9] = rho;
844     vals[10] = vertices;
845     vals[11] = EtaSeed - scEta;
846     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
847     vals[13] = ESeed/SCRawEnergy;
848     vals[14] = E3x3Seed/ESeed;
849     vals[15] = E5x5Seed/ESeed;
850     vals[16] = see;
851     vals[17] = spp;
852     vals[18] = sep;
853     vals[19] = EMaxSeed/ESeed;
854     vals[20] = E2ndSeed/ESeed;
855     vals[21] = ETopSeed/ESeed;
856     vals[22] = EBottomSeed/ESeed;
857     vals[23] = ELeftSeed/ESeed;
858     vals[24] = ERightSeed/ESeed;
859     vals[25] = E2x5MaxSeed/ESeed;
860     vals[26] = E2x5TopSeed/ESeed;
861     vals[27] = E2x5BottomSeed/ESeed;
862     vals[28] = E2x5LeftSeed/ESeed;
863     vals[29] = E2x5RightSeed/ESeed;
864     vals[30] = GsfTrackPIn;
865     vals[31] = fbrem;
866     vals[32] = Charge;
867     vals[33] = EoP;
868     vals[34] = TrackMomentumError/GsfTrackPIn;
869     vals[35] = detaIn;
870     vals[36] = dphiIn;
871     vals[37] = detaCalo;
872     vals[38] = dphiCalo;
873     vals[39] = GsfTrackChiSqr;
874     vals[40] = KFTrackNLayers;
875     vals[41] = ElectronEnergyOverPout;
876     vals[42] = PreShowerOverRaw;
877     }
878    
879     // Now evaluating the regression
880     double regressionResult = 0;
881    
882     if (fVersionType == kWithTrkVar) {
883     if (fabs(scEta) <= 1.479) regressionResult = SCRawEnergy * forestCorrection_eb->GetResponse(vals);
884     else regressionResult = (SCRawEnergy*(1+PreShowerOverRaw)) * forestCorrection_ee->GetResponse(vals);
885     }
886    
887    
888     //print debug
889     if (printDebug) {
890     if (scEta <= 1.479) {
891     std::cout << "Barrel :";
892     for (uint v=0; v < 50; ++v) std::cout << vals[v] << ", ";
893     std::cout << "\n";
894     }
895     else {
896     std::cout << "Endcap :";
897     for (uint v=0; v < 43; ++v) std::cout << vals[v] << ", ";
898     std::cout << "\n";
899     }
900     std::cout << "SCRawEnergy = " << SCRawEnergy << " : PreShowerOverRaw = " << PreShowerOverRaw << std::endl;
901     std::cout << "regression energy = " << regressionResult << std::endl;
902     }
903    
904     // Cleaning up and returning
905     delete[] vals;
906     return regressionResult;
907     }
908    
909    
910    
911    
912     double ElectronEnergyRegression::regressionUncertaintyWithTrkVarV2( std::vector<double> &inputvars,
913     bool printDebug)
914     {
915     // Checking if instance has been initialized
916     if (fIsInitialized == kFALSE) {
917     printf("ElectronEnergyRegression instance not initialized !!!");
918     return 0;
919     }
920    
921     //Check that inputvars vector has the right number of inputs
922     assert(inputvars.size() == 46);
923    
924     // Checking if fVersionType is correct
925     assert(fVersionType == kWithTrkVar);
926    
927     //assign variables from inputvars array to named variables
928     double SCRawEnergy = inputvars[0];
929     double scEta = inputvars[1];
930     double scPhi = inputvars[2];
931     double R9 = inputvars[3];
932     double etawidth = inputvars[4];
933     double phiwidth = inputvars[5];
934     double NClusters = inputvars[6];
935     double HoE = inputvars[7];
936     double rho = inputvars[8];
937     double vertices = inputvars[9];
938     double EtaSeed = inputvars[10];
939     double PhiSeed = inputvars[11];
940     double ESeed = inputvars[12];
941     double E3x3Seed = inputvars[13];
942     double E5x5Seed = inputvars[14];
943     double see = inputvars[15];
944     double spp = inputvars[16];
945     double sep = inputvars[17];
946     double EMaxSeed = inputvars[18];
947     double E2ndSeed = inputvars[19];
948     double ETopSeed = inputvars[20];
949     double EBottomSeed = inputvars[21];
950     double ELeftSeed = inputvars[22];
951     double ERightSeed = inputvars[23];
952     double E2x5MaxSeed = inputvars[24];
953     double E2x5TopSeed = inputvars[25];
954     double E2x5BottomSeed = inputvars[26];
955     double E2x5LeftSeed = inputvars[27];
956     double E2x5RightSeed = inputvars[28];
957     double IEtaSeed = inputvars[29];
958     double IPhiSeed = inputvars[30];
959     double EtaCrySeed = inputvars[31];
960     double PhiCrySeed = inputvars[32];
961     double PreShowerOverRaw = inputvars[33];
962     double GsfTrackPIn = inputvars[34];
963     double fbrem = inputvars[35];
964     double Charge = inputvars[36];
965     double EoP = inputvars[37];
966     double TrackMomentumError = inputvars[38];
967     double detaIn = inputvars[39];
968     double dphiIn = inputvars[40];
969     double detaCalo = inputvars[41];
970     double dphiCalo = inputvars[42];
971     double GsfTrackChiSqr = inputvars[43];
972     double KFTrackNLayers = inputvars[44];
973     double ElectronEnergyOverPout = inputvars[45];
974    
975     float *vals = (fabs(scEta) <= 1.479) ? new float[50] : new float[43];
976     if (fabs(scEta) <= 1.479) { // Barrel
977     vals[0] = SCRawEnergy;
978     vals[1] = scEta;
979     vals[2] = scPhi;
980     vals[3] = R9;
981     vals[4] = E5x5Seed/SCRawEnergy;
982     vals[5] = etawidth;
983     vals[6] = phiwidth;
984     vals[7] = NClusters;
985     vals[8] = HoE;
986     vals[9] = rho;
987     vals[10] = vertices;
988     vals[11] = EtaSeed - scEta;
989     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
990     vals[13] = ESeed/SCRawEnergy;
991     vals[14] = E3x3Seed/ESeed;
992     vals[15] = E5x5Seed/ESeed;
993     vals[16] = see;
994     vals[17] = spp;
995     vals[18] = sep;
996     vals[19] = EMaxSeed/ESeed;
997     vals[20] = E2ndSeed/ESeed;
998     vals[21] = ETopSeed/ESeed;
999     vals[22] = EBottomSeed/ESeed;
1000     vals[23] = ELeftSeed/ESeed;
1001     vals[24] = ERightSeed/ESeed;
1002     vals[25] = E2x5MaxSeed/ESeed;
1003     vals[26] = E2x5TopSeed/ESeed;
1004     vals[27] = E2x5BottomSeed/ESeed;
1005     vals[28] = E2x5LeftSeed/ESeed;
1006     vals[29] = E2x5RightSeed/ESeed;
1007     vals[30] = GsfTrackPIn;
1008     vals[31] = fbrem;
1009     vals[32] = Charge;
1010     vals[33] = EoP;
1011     vals[34] = TrackMomentumError/GsfTrackPIn;
1012     vals[35] = detaIn;
1013     vals[36] = dphiIn;
1014     vals[37] = detaCalo;
1015     vals[38] = dphiCalo;
1016     vals[39] = GsfTrackChiSqr;
1017     vals[40] = KFTrackNLayers;
1018     vals[41] = ElectronEnergyOverPout;
1019     vals[42] = IEtaSeed;
1020     vals[43] = IPhiSeed;
1021     vals[44] = ((int) IEtaSeed)%5;
1022     vals[45] = ((int) IPhiSeed)%2;
1023     vals[46] = (abs(IEtaSeed)<=25)*(((int)IEtaSeed)%25) + (abs(IEtaSeed)>25)*(((int) (IEtaSeed-25*abs(IEtaSeed)/IEtaSeed))%20);
1024     vals[47] = ((int) IPhiSeed)%20;
1025     vals[48] = EtaCrySeed;
1026     vals[49] = PhiCrySeed;
1027     }
1028    
1029     else { // Endcap
1030     vals[0] = SCRawEnergy;
1031     vals[1] = scEta;
1032     vals[2] = scPhi;
1033     vals[3] = R9;
1034     vals[4] = E5x5Seed/SCRawEnergy;
1035     vals[5] = etawidth;
1036     vals[6] = phiwidth;
1037     vals[7] = NClusters;
1038     vals[8] = HoE;
1039     vals[9] = rho;
1040     vals[10] = vertices;
1041     vals[11] = EtaSeed - scEta;
1042     vals[12] = atan2(sin(PhiSeed-scPhi),cos(PhiSeed-scPhi));
1043     vals[13] = ESeed/SCRawEnergy;
1044     vals[14] = E3x3Seed/ESeed;
1045     vals[15] = E5x5Seed/ESeed;
1046     vals[16] = see;
1047     vals[17] = spp;
1048     vals[18] = sep;
1049     vals[19] = EMaxSeed/ESeed;
1050     vals[20] = E2ndSeed/ESeed;
1051     vals[21] = ETopSeed/ESeed;
1052     vals[22] = EBottomSeed/ESeed;
1053     vals[23] = ELeftSeed/ESeed;
1054     vals[24] = ERightSeed/ESeed;
1055     vals[25] = E2x5MaxSeed/ESeed;
1056     vals[26] = E2x5TopSeed/ESeed;
1057     vals[27] = E2x5BottomSeed/ESeed;
1058     vals[28] = E2x5LeftSeed/ESeed;
1059     vals[29] = E2x5RightSeed/ESeed;
1060     vals[30] = GsfTrackPIn;
1061     vals[31] = fbrem;
1062     vals[32] = Charge;
1063     vals[33] = EoP;
1064     vals[34] = TrackMomentumError/GsfTrackPIn;
1065     vals[35] = detaIn;
1066     vals[36] = dphiIn;
1067     vals[37] = detaCalo;
1068     vals[38] = dphiCalo;
1069     vals[39] = GsfTrackChiSqr;
1070     vals[40] = KFTrackNLayers;
1071     vals[41] = ElectronEnergyOverPout;
1072     vals[42] = PreShowerOverRaw;
1073     }
1074    
1075     // Now evaluating the regression
1076     double regressionResult = 0;
1077    
1078     if (fVersionType == kWithTrkVar) {
1079     if (fabs(scEta) <= 1.479) regressionResult = SCRawEnergy * forestUncertainty_eb->GetResponse(vals);
1080     else regressionResult = (SCRawEnergy*(1+PreShowerOverRaw)) * forestUncertainty_ee->GetResponse(vals);
1081     }
1082    
1083     //print debug
1084     if (printDebug) {
1085     if (scEta <= 1.479) {
1086     std::cout << "Barrel :";
1087     for (uint v=0; v < 50; ++v) std::cout << vals[v] << ", ";
1088     std::cout << "\n";
1089     }
1090     else {
1091     std::cout << "Endcap :";
1092     for (uint v=0; v < 43; ++v) std::cout << vals[v] << ", ";
1093     std::cout << "\n";
1094     }
1095     std::cout << "SCRawEnergy = " << SCRawEnergy << " : PreShowerOverRaw = " << PreShowerOverRaw << std::endl;
1096     std::cout << "regression energy uncertainty = " << regressionResult << std::endl;
1097     }
1098    
1099    
1100     // Cleaning up and returning
1101     delete[] vals;
1102     return regressionResult;
1103     }
1104