CalorimeterEnergyCorrelatorCommonTools/EventShapeQA/EventShapeQA.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef EVENTSHAPEQA_H
0004 #define EVENTSHAPEQA_H
0005
0006 //fun4all
0007 #include <fun4all/SubsysReco.h>
0008 #include <fun4all/Fun4AllBase.h>
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010
0011 //phool
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/getClass.h>
0014 #include <phool/PHNodeOperation.h>
0015
0016 //Calo towers
0017 #include <calobase/TowerInfoContainer.h>
0018 #include <calobase/TowerInfoContainerv1.h>
0019 #include <calobase/TowerInfoContainerv2.h>
0020 #include <calobase/TowerInfov2.h>
0021 #include <calobase/TowerInfov1.h>
0022 #include <calobase/TowerInfo.h>
0023 #include <calobase/RawTowerDefs.h>
0024 #include <calobase/RawTowerContainer.h>
0025 #include <calobase/RawTowerGeomContainer.h>
0026 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0027
0028 //G4 objects
0029 #include <g4main/PHG4Particle.h>
0030 #include <g4main/PHG4Hit.h>
0031 #include <g4main/PHG4TruthInfoContainer.h>
0032
0033 //jetbase objetcts
0034 #include <jetbase/JetContainer.h>
0035 #include <jetbase/JetContainerv1.h>
0036 #include <jetbase/JetMapv1.h>
0037 #include <jetbase/JetMap.h>
0038 #include <jetbase/Jetv1.h>
0039 #include <jetbase/Jetv2.h>
0040 #include <jetbase/Jet.h>
0041 #include <jetbase/TowerJetInput.h>
0042 #include <jetbase/ClusterJetInput.h>
0043
0044 //fastjet
0045 #include <fastjet/PseudoJet.hh>
0046 #include <fastjet/ClusterSequence.hh>
0047
0048 //vertex stuff
0049 #include <globalvertex/GlobalVertex.h>
0050 #include <globalvertex/GlobalVertexMap.h>
0051
0052 //c++
0053 #include <map>
0054 #include <utility>
0055 #include <vector>
0056 #include <string>
0057 #include <math.h>
0058 #include <array>
0059 #include <iostream>
0060
0061 //root
0062 #include <TH1.h>
0063 #include <TH2.h>
0064 #include <TFile.h>
0065 #include <TTree.h>
0066 #include <TInterpreter.h>
0067
0068 #define PI 3.14159265358979
0069 class PHCompositeNode;
0070 enum CompDataTypes
0071 {
0072     notComp = -1,
0073     G4TruthMap = 0,
0074     G4TruthPrimaryParticleRange = 1,
0075     PHHepMCGenEvent = 2,
0076     AllCalTower = 3,
0077     EMCALTower = 4,
0078     IHCALTower = 5,
0079     OHCALTower = 6,
0080     TopoCluster = 7,
0081     AllCALTowerCluster  = 8,
0082     EMCALTowerCluster = 9,
0083     IHCALTowerCluster = 10,
0084     OHCALTowerCluster = 11
0085 };
0086 enum JetDataTypes
0087 {
0088     notJet = -1,
0089     G4TruthJets = 0,
0090     TowerJets = 1,
0091     ClusterJets = 2
0092 };
0093
0094 class HitPlots
0095 {
0096   public:
0097     HitPlots(std::string input_type="G4 Particle Truth", bool jet=false, bool shifted=false)
0098     {
0099         isJet=jet;
0100         isShifted=shifted;
0101         inputType=input_type;
0102         convertToEnum(inputType, isJet);
0103         convertToLabel(this->ComponentType, this->JetType, isJet);
0104         makePlots();
0105     }
0106     void convertToEnum(std::string it, bool iJ)
0107     {
0108         if(iJ)
0109         {
0110             if(it.find("G4") != std::string::npos || it.find("Truth") != std::string::npos) this->JetType = JetDataTypes::G4TruthJets;
0111             else if(it.find("ower" ) != std::string::npos) this->JetType = JetDataTypes::TowerJets;
0112             else if(it.find("luster") != std::string::npos) this->JetType = JetDataTypes::ClusterJets;
0113             else this->JetType = JetDataTypes::notJet;
0114             this->ComponentType = CompDataTypes::notComp;
0115         }
0116         else{
0117             if(it.find("G4") != std::string::npos && it.find("Primary") == std::string::npos) this->ComponentType = CompDataTypes::G4TruthMap;
0118             else if( it.find("G4") != std::string::npos) this->ComponentType = CompDataTypes::G4TruthPrimaryParticleRange;
0119             else if ( it.find("HepMC") != std::string::npos) this->ComponentType = CompDataTypes::PHHepMCGenEvent;
0120             else if ( it.find("ower") != std::string::npos)
0121             {
0122                 if( it.find("luster") == std::string::npos)
0123                 {
0124                     if(       it.find("All")!= std::string::npos || it.find("ALL" != std::string::npos) this->ComponentType = CompDataTypes::AllCalTower;
0125                     else if ( it.find("EM") != std::string::npos || it.find("em") != std::string::npos) this->ComponentType = CompDataTypes::EMCALTower;
0126                     else if ( it.find("IH") != std::string::npos || it.find("ih") != std::string::npos) this->ComponentType = CompDataTypes::IHCALTower;
0127                     else if ( it.find("OH") != std::string::npos || it.find("oh") != std::string::npos) this->ComponentType = CompDataTypes::OHCALTower;
0128                 }
0129                 else
0130                 {
0131                     if(       it.find("All")!= std::string::npos || it.find("ALL" != std::string::npos) this->ComponentType = CompDataTypes::AllCalTower;
0132                     else if ( it.find("EM") != std::string::npos || it.find("em") != std::string::npos) this->ComponentType = CompDataTypes::EMCALTower;
0133                     else if ( it.find("IH") != std::string::npos || it.find("ih") != std::string::npos) this->ComponentType = CompDataTypes::IHCALTower;
0134                     else if ( it.find("OH") != std::string::npos || it.find("oh") != std::string::npos) this->ComponentType = CompDataTypes::OHCALTower;
0135
0136                 }
0137             }
0138             else if ( it.find("luster") != std::string::npos ) this->ComponentType = CompDataTypes::TopoCluster;
0139             else this->ComponentType = CompDataTypes::notComp;
0140             this->JetType = JetDataTypes::nonJet;
0141             }
0142             return;
0143     }
0144
0145     void convertToLabel(CompDataTypes ct, JetDataTypes jt, bool iJ=false)
0146     {
0147         if(!iJ)
0148         {
0149             switch (ct)
0150             {
0151                 case CompDataTypes::G4TruthMap :
0152                     this->generic_name_label+="G4truth_map_";
0153                     this->generic_Title_label+="over G4 Truth Map Components ";
0154                     break;
0155                 case CompDataTypes::G4TruthPrimaryParticleRange :
0156                     this->generic_name_label+="G4Truth_PrimaryParticle_Range_;
0157                     this->generic_Title_label+="over G4 Truth Primary Particle Components ;
0158                     break;
0159                 case CompDataTypes::PHHepMCGenEvent  :
0160                     this->generic_name_label+="HepMC_map_";
0161                     this->generic_Title_label+="over HepMC map ";
0162                     break;
0163                 case CompDataTypes:: AllCALTower  :
0164                     this->generic_name_label+="AllCAL_Tower_";
0165                     this->generic_Title_label+="over Combined Calorimeter Towers ";
0166                     break;
0167                 case CompDataTypes:: EMCALTower  :
0168                     this->generic_name_label+="EMCAL_Tower_";
0169                     this->generic_Title_label+="over recombined EMCAL Towers ";
0170                     break;
0171                 case CompDataTypes:: IHCALTower  :
0172                     this->generic_name_label+="IHCAL_Tower_";
0173                     this->generic_Title_label+="over IHCAL Towers ";
0174                     break;
0175                 case CompDataTypes:: OHCALTower :
0176                     this->generic_name_label+="OHCAL_Tower_";
0177                     this->generic_Title_label+="over OHCAL Towers ";
0178                     break;
0179                 case CompDataTypes::TopoCluster:
0180                     this->generic_name_label+="Cluster_";
0181                     this->generic_Title_label+="over Clusters ";
0182                     break;
0183                 case CompDataTypes::ALLCALTowerCluster :
0184                     this->generic_name_label+="CAL_Cluster_Towers_";
0185                     this->generic_Title_label+="over combined Calorimeter Towers identified from clusters";
0186                     break;
0187                 case CompDataTypes::EMCALTowerCluster :
0188                     this->generic_name_label+="EMCAL_Cluster_Towers_";
0189                     this->generic_Title_label+="over retowered EMCAL Towers identified from clusters";
0190                     break;
0191                 case CompDataTypes::IHCALTowerCluster :
0192                     this->generic_name_label+="IHCAL_Cluster_Towers_";
0193                     this->generic_Title_label+="over IHCAL Towers identified from clusters";
0194                     break;
0195                 case CompDataTypes::OHCALTowerCluster :
0196                     this->generic_name_label+="OHCAL_Cluster_Towers_";
0197                     this->generic_Title_label+="over OHCAL Towers identified from clusters";
0198                     break;
0199
0200             }
0201         }
0202         else
0203         {
0204             switch (jt){
0205                 case JetDataTypes::G4TruthJets :
0206                     this->generic_name_label+="G4truth_jets_";
0207                     this->generic_Title_label+="over G4 Truth Jets ";
0208                     break;
0209                 case JetDataTypes::TowerJets :
0210                     this->generic_name_label+="Tower_jets_";
0211                     this->generic_Title_label+="over Tower Jets ";
0212                     break;
0213                 case JetDataTypes::ClusterJets :
0214                     this->generic_name_label+="Cluster_jets_";
0215                     this->generic_Title_label+="over Cluster Jets "
0216                     break;
0217             }
0218         }
0219         if(this->isShifted){
0220             this->generic_name_label+="relative_axis_";
0221             this->generic_Title_label+="relative to Lead axis ";
0222         }
0223         else
0224         {
0225             this->generic_name_label+="zero_axis_";
0226             this->generic_Title_label+="relative to detector axis ";
0227         }
0228         return;
0229     }
0230     void makePlots()
0231     {
0232         float pm=-PI, px=PI;
0233         float em=-1.2, ex=1.2;
0234         float Em=-0.5, Ex=99.5;
0235         float Rm=0, Rx=4.5;
0236         float xjm=0., xjx=1.;
0237         float Aj=0., Ajx=1.;
0238         bool istruth=true;
0239         xj=new TH1F(Form("%s_xj", this->generic_name_label.c_str()), Form("x_{j} %s ; x_{j}; Counts ", this->generic_Title_label.c_str()), 100, xjm, xjx);
0240         Aj=new TH1F(Form("%s_Aj", this->generic_name_label.c_str()), Form("A_{j} %s ; A_{j}; Counts ", this->generic_Title_label.c_str()), 100, Ajm, Ajx);
0241         E=new TH1F(Form("%s_E", this->generic_name_label.c_str()), Form("E %s ; E [GeV]; Counts ", this->generic_Title_label.c_str()), 100, Em, Ex);
0242         ET=new TH1F(Form("%s_ET", this->generic_name_label.c_str()), Form("E_{T} %s ; E_{T} [GeV]; Counts ", this->generic_Title_label.c_str()), 100, Em, Ex);
0243         R=new TH1F(Form("%s_R", this->generic_name_label.c_str()), Form("#Delta R from central axis %s ; #Delta R; Counts ", this->generic_Title_label.c_str()), 100, Rm, Rx);
0244         phi=new TH1F(Form("%s_phi", this->generic_name_label.c_str()), Form("varphi %s ; #varphi; Counts ", this->generic_Title_label.c_str()), 100, pm, px);
0245         eta=new TH1F(Form("%s_eta", this->generic_name_label.c_str()), Form("#eta %s ; #eta; Counts ", this->generic_Title_label.c_str()), 100, em, ex);
0246             eta_phi_hit = new TH2F(Form("%s_phi_eta_hit", this->generic_name_label.c_str()), Form("Hit map %s; #eta; #varphi", this->generic_Title_label.c_str()), 100, em, ex, 100, pm, px);
0247             eta_phi_E = new TH2F(Form("%s_phi_eta_E", this->generic_name_label.c_str()), Form("Energy Deposition %s; #eta; #varphi", this->generic_Title_label.c_str()), 100, em, ex, 100, pm, px);
0248             E_R_hit = new TH2F(Form("%s_E_R_hit", this->generic_name_label.c_str()), Form("Hit map %s;E [GeV]; #Delta R", this->generic_Title_label.c_str()), 100, Em, Ex, 100, Rm, Rx);
0249         return;
0250     }
0251     TH1F* xj, *Aj, *phi, *eta, *E, *ET, *R;
0252     TH2F* eta_phi_hit, *eta_phi_E, *E_R_hit;
0253     CompDataTypes ComponentType = CompDataTypes::notComp;
0254     JetDataTypes JetType = JetDataTypes::notJet;
0255     bool isJet = false;
0256     bool isShifted=false;
0257     std::string generic_name_label="h_";
0258     std::string generic_Title_label="";
0259 };
0260
0261
0262 class EventShapeQA : public SubsysReco
0263 {
0264  public:
0265
0266     EventShapeQA(int verb = 0, const std::string &name = "EventShapeQA");
0267
0268     ~EventShapeQA() override;
0269
0270     int Init(PHCompositeNode *topNode) override;
0271
0272     int InitRun(PHCompositeNode *topNode) override {return Fun4AllReturnCodes::EVENT_OK;}
0273
0274     int process_event(PHCompositeNode *topNode) override;
0275
0276     int ResetEvent(PHCompositeNode *topNode) override {return Fun4AllReturnCodes::EVENT_OK;}
0277
0278     int EndRun(const int runnumber) override {return Fun4AllReturnCodes::EVENT_OK;}
0279
0280     int End(PHCompositeNode *topNode) override;
0281
0282     int Reset(PHCompositeNode * /*topNode*/) override {return Fun4AllReturnCodes::EVENT_OK;}
0283
0284     void Print(const std::string &what = "ALL") const override;
0285     float PhiWrap(float phi)
0286     {
0287         if(phi > PI) phi = 2*PI - phi;
0288         else if (phi < - PI ) phi = - 2 *PI - phi;
0289         return phi;
0290     }
0291     float calcR(float phi1, float phi2, float eta1, float eta2)
0292     {
0293         float dphi = phi1-phi2;
0294         dphi = PhiWrap(dphi);
0295         float deta = eta1 - eta2;
0296         float R = std::sqrt(std::pow(dphi, 2) + std::pow(deta, 2));
0297         return R;
0298     }
0299
0300  private:
0301     void doPHG4Analysis(std::pair<std::pair<CompDataTypes, JetDataTypes>, PHCompositeNode*>>)
0302     void doPHG4JetAnalysis(std::pair<std::pair<CompDataTypes, JetDataTypes>, PHCompositeNode*>>)
0303     void doHepMCAnalysis(std::pair<std::pair<CompDataTypes, JetDataTypes>, PHCompositeNode*>>)
0304     void doCaloAnalysis(std::pair<std::pair<CompDataTypes, JetDataTypes>, PHCompositeNode*>>)
0305     void doCaloJetAnalysis(std::pair<std::pair<CompDataTypes, JetDataTypes>, PHCompositeNode*>>)
0306     int verbosity = 0, n_evt = 0;
0307     //particles in reference to the proper calorimeter axis
0308     std::vector<std::pair<std::pair<CompDataTypes, JetDataTypes> , HitPlots*>>* calo_jet_zero_axis=NULL;
0309     //jets in reference to the proper calorimeter axis
0310     std::vector<std::pair<std::pair<CompDataTypes, JetDataTypes> , HitPlots*>>* calo_jet_shift_axis=NULL;
0311     std::vector<std::pair<std::pair<CompDataTypes, JetDataTypes> ,PHCompositeNode*>> nodes {};
0312
0313 };
0314
0315 #endif // EVENTSHAPEQA_H