LargeRLCaloENC/src/HelperStructs.h

0001 #ifndef __HELPERSTRUCTS_H__
0002 #define __HELPERSTRUCTS_H__
0003 #include <TH1.h>
0004 #include <TH2.h>
0005
0006 #include <vector>
0007 #include <string>
0008 #include <map>
0009 #include <math.h>
0010
0011 class TowerOutput
0012 {
0013     public:
0014         TowerOutput(float tr=0)
0015         {
0016             this->threshold=tr;
0017         }
0018         ~TowerOutput(){};
0019         void AddE2CValues(float e2c, float rl){
0020             int index=-1;
0021             //see if we already have an entry in this vertex
0022             for(int i = 0; i<(int) RL.size(); i++){
0023                 if(rl == RL.at(i) ){
0024                     index = i;
0025                     break;
0026                 }
0027                 else continue;
0028             }
0029             if(index > -1 ){
0030                 //just add to the value if we already have other elements at this fixed value
0031                 this->E2C.at(index)+=e2c;
0032             }
0033             else{
0034                 this->E2C.push_back(e2c);
0035                 this->RL.push_back(rl);
0036             }
0037             return;
0038         }
0039         void AddE3CValues(float e3c, std::array<float, 3> dist){
0040             int index=-1;
0041             if(RL_RM_RS.size() > 0 ){
0042             for(int i = 0; i<(int)RL_RM_RS.size(); i++){
0043                 if(i >= (int) RL_RM_RS.size()) break;
0044                 if(dist == RL_RM_RS.at(i)){
0045                     index = i;
0046                     break;
0047                 }
0048                 else continue;
0049             }}
0050             if(index > -1 ){
0051                 E3C_full_shape.at(index)+=e3c;
0052             }
0053             else{
0054                 E3C_full_shape.push_back(e3c);
0055                 RL_RM_RS.push_back(dist);
0056             }
0057             return;
0058         }
0059         void CalculateFlatE3C()
0060         {
0061             std::map<float,float> e3c;
0062             for(int i= 0; i<(int)RL_RM_RS.size(); i++){
0063                 float rl=RL_RM_RS.at(i)[0];
0064                 if(e3c.find(rl) != e3c.end()){
0065                     e3c[rl]+=E3C_full_shape.at(i);
0066                 }
0067                 else{
0068                     e3c[rl]=E3C_full_shape.at(i);
0069                 }
0070             }
0071             for(auto r:RL)
0072             {
0073                 if(e3c.find(r) != e3c.end()) E3C.push_back(e3c[r]);
0074             } //just keeps the ordering fixed
0075             return;
0076         }
0077         bool Merge(std::vector<TowerOutput*> tos)
0078         {
0079             std::map<float, float> e2c;
0080             std::map<std::array<float, 3>, float> e3c;
0081             for(auto to:tos){
0082                 if(this->threshold != to->threshold){
0083                     return false;
0084                 }
0085                 for(int i = 0; i < (int)to->RL.size(); i++){
0086                     if(e2c.find(to->RL.at(i)) !=e2c.end())
0087                     {
0088                         e2c[to->RL.at(i)]+=to->E2C.at(i);
0089                     }
0090                     else e2c[to->RL.at(i)]=to->E2C.at(i);
0091                 }
0092                 for(int i = 0; i < (int)to->RL_RM_RS.size(); i++){
0093                     if(e3c.find(to->RL_RM_RS.at(i)) != e3c.end())
0094                     {
0095                         e3c[to->RL_RM_RS.at(i)]+=to->E3C_full_shape.at(i);
0096                     }
0097                     else  e3c[to->RL_RM_RS.at(i)]+=to->E3C_full_shape.at(i);
0098                 }
0099             }
0100
0101                 //create a merge function that preserves indexing of unique values
0102             for(auto e:e2c)
0103             {
0104                 this->AddE2CValues(e.second, e.first);
0105             }
0106             for(auto e:e3c)
0107             {
0108                     this->AddE3CValues( e.second, e.first);
0109                 }
0110             //  this->CalculateFlatE3C();
0111             return true;
0112         }
0113
0114
0115         void Normalize(float energy)
0116         {
0117             for(auto e2c: this->E2C)
0118             {
0119                 e2c=e2c*(float)std::pow(energy,-2);
0120             }
0121             for(auto e3c: this->E3C)
0122             {
0123                 e3c=e3c*(float)std::pow(energy, -3);
0124             }
0125             for(auto e3c: this->E3C_full_shape)
0126             {
0127                 e3c=e3c*(float)std::pow(energy, -3);
0128             }
0129             return;
0130         }
0131         float threshold = 0.;
0132         std::vector<float> RL;
0133         std::vector<std::array<float, 3>> RL_RM_RS;
0134         std::vector<float> E2C;
0135         std::vector<float> E3C;
0136         std::vector<float> E3C_full_shape;
0137
0138
0139 };
0140
0141 class StrippedDownTower
0142 {
0143     public:
0144         enum Regions{
0145             Full,
0146             Towards,
0147             Away,
0148             TransverseMax,
0149             TransverseMin
0150         };
0151         enum Calorimeter{
0152             All,
0153             EMCAL,
0154             IHCAL,
0155             OHCAL,
0156             TRUTH
0157         };
0158         StrippedDownTower( float threshold, Calorimeter which_calo=Calorimeter::All, Regions region_tag = Regions::Full){
0159             //initialize the stripped down tower object
0160             this->tag = region_tag;
0161             this->cal = which_calo;
0162             this->RegionOutput=new TowerOutput(threshold);
0163             this->FullOutput=new TowerOutput(threshold);
0164
0165         }
0166         ~StrippedDownTower(){};
0167         float r     =0.;
0168         float phi   =0.;
0169         float eta   =0.;
0170         float E     =0.;
0171         Regions tag =Regions::Full;
0172         Calorimeter cal =Calorimeter::All;
0173         TowerOutput* RegionOutput;
0174         TowerOutput* FullOutput;
0175 };
0176 class StrippedDownTowerWithThresholds
0177 {
0178     public:
0179         enum Regions{
0180             Full,
0181             Towards,
0182             Away,
0183             TransverseMax,
0184             TransverseMin
0185         };
0186         enum Calorimeter{
0187             All,
0188             EMCAL,
0189             IHCAL,
0190             OHCAL
0191         };
0192         StrippedDownTowerWithThresholds( std::vector<float> thresholds, Calorimeter which_calo=Calorimeter::All, Regions region_tag = Regions::Full){
0193             //initialize the stripped down tower object
0194             this->tag = region_tag;
0195             this->cal = which_calo;
0196             this->RegionOutput=new std::vector<TowerOutput*>;
0197             this->FullOutput=new std::vector<TowerOutput*>;
0198             for(float threshold:thresholds){
0199                 TowerOutput* to=new TowerOutput(threshold);
0200                 TowerOutput* tf=new TowerOutput(threshold);
0201                 this->RegionOutput->push_back(to);
0202                 this->FullOutput->push_back(tf);
0203             }
0204
0205         }
0206         ~StrippedDownTowerWithThresholds(){};
0207         int getThresholdIndex(float threshold, bool RegionOrFull)
0208         {
0209             int index=-1;
0210             if(RegionOrFull){
0211                 for(int i=0; i< (int)this->FullOutput->size(); i++){
0212                     float thresh_temp=this->FullOutput->at(i)->threshold;
0213                     if(thresh_temp == threshold){
0214                         index = i;
0215                         break;
0216                     }
0217                     else continue;
0218                 }
0219             }
0220             else if(!RegionOrFull){
0221                 for(int i=0; i< (int)this->RegionOutput->size(); i++){
0222                     float thresh_temp=this->RegionOutput->at(i)->threshold;
0223                     if(thresh_temp == threshold){
0224                         index = i;
0225                         break;
0226                     }
0227                     else continue;
0228                 }
0229             }
0230             return index;
0231         }
0232
0233
0234
0235         float r     =0.;
0236         float phi   =0.;
0237         float eta   =0.;
0238         float E     =0.;
0239         Regions tag =Regions::Full;
0240         Calorimeter cal =Calorimeter::All;
0241         std::vector<TowerOutput*>* RegionOutput;
0242         std::vector<TowerOutput*>* FullOutput;
0243
0244 };
0245 struct DijetQATypePlots{
0246     DijetQATypePlots(){
0247         bad_occ_em_oh_rat=new TH2F("h_bad_occupancy_EM_OH_rat", "Occupancy for events that fail the OHCAL ratio cut; #% Towers #geq 10 MeV EMCAL; #% Towers > 10 MeV OHCAL; N_{Evts}", 100, -0.005, 0.995, 100, -0.005, 0.995);
0248         bad_occ_em_h_rat=new TH2F("h_bad_occupancy_EM_H_rat", "Occupancy for events that fail OHCAL ratio cut; #% Towers #geq 10 MeV EMCAL; (#% Towers #geq 10 MeV OHCAL + #% Towers #geq 10 MeV IHCAL)/2", 100, -0.005, 0.995, 100, -0.005, 0.995);
0249         bad_occ_em_oh=new TH2F("h_bad_occupancy_EM_OH", "Occupancy for events that pass  OHCAL ratio cut but otherwise fail; #% Towers #geq 10 MeV EMCAL; #% Towers > 10 MeV OHCAL; N_{Evts}", 100, -0.005, 0.995, 100, -0.005, 0.995);
0250         bad_occ_em_h=new TH2F("h_bad_occupancy_EM_H", "Occupancy for events that pass OHCAL ratio cut but otherwise fail; #% Towers #geq 10 MeV EMCAL; (#% Towers #geq 10 MeV OHCAL + #% Towers #geq 10 MeV IHCAL)/2", 100, -0.005, 0.995, 100, -0.005, 0.995);
0251         ohcal_bad_hits=new TH2F("h_ohcal_bad_hits", "Energy depositon for events that fail the OHCAL ratio cut; #eta; #varphi; E [GeV]", 24, -1.1, 1.1, 64, -PI, PI);
0252         emcal_occup=new TH1F("h_EMCAL_occupancy", "Occupancy of EMCAL all events; #% Towers #geq 10 MeV EMCAL; N_{evts}", 100, -0.005, 0.995);
0253         ihcal_occup=new TH1F("h_IHCAL_occupancy", "Occupancy of IHCAL all events; #% Towers #geq 10 MeV IHCAL; N_{evts}", 100, -0.005, 0.995);
0254         ohcal_occup=new TH1F("h_OHCAL_occupancy", "Occupancy of OHCAL all events; #% Towers #geq 10 MeV OHCAL; N_{evts}", 100, -0.005, 0.995);
0255         ohcal_rat_occup=new TH2F("h_OHCAL_rat", "Occupancy of OHCAL as function of OHCAL ratio; #frac{E_{OHCAL}}{E_{ALL CALS}}; #% Towers #geq 500 MeV; N_{evts}", 100, -0.005, 0.995, 100, -0.005, 0.995);
0256     }
0257     TH2F* bad_occ_em_oh_rat=nullptr;
0258     TH2F* bad_occ_em_h_rat=nullptr;
0259     TH2F* bad_occ_em_oh=nullptr;
0260     TH2F* bad_occ_em_h=nullptr;
0261     TH2F* ohcal_bad_hits=nullptr;
0262     TH2F* ohcal_rat_occup=nullptr;
0263     TH1F* emcal_occup=nullptr;
0264     TH1F* ihcal_occup=nullptr;
0265     TH1F* ohcal_occup=nullptr;
0266     TH1F* ohcal_rat_h=nullptr;
0267     std::vector<TH2F*> QA_2D_plots {bad_occ_em_oh_rat, bad_occ_em_h_rat, bad_occ_em_oh, bad_occ_em_h, ohcal_bad_hits, ohcal_rat_occup};
0268     std::vector<TH1F*> QA_1D_plots {emcal_occup, ihcal_occup, ohcal_occup};
0269     };
0270 #endif