---

EtMissTest.cpp

Go to the documentation of this file.

#include <Utilities/Configuration/interface/Architecture.h>

#include "CARF/G3Event/interface/G3EventProxy.h"
#include "Utilities/Notification/interface/Observer.h"

// Calorimetry stuff
#include "Calorimetry/CaloRecHit/interface/CaloRecHit.h"
#include "Calorimetry/EcalPlusHcalTower/interface/EcalPlusHcalTower.h"
#include "Calorimetry/EcalPlusHcalTower/src/EcalPlusHcalTowerBuilder.h"
#include "Calorimetry/CaloCommon/interface/CaloItr.h"
#include "Calorimetry/CaloDetector/interface/CellID.h"

// reconstructor stuff
#include "CARF/Reco/interface/RecItr.h"
#include "CARF/Reco/interface/Reconstructor.h"

#include "Utilities/Notification/interface/PackageInitializer.h"
#include <Utilities/UI/interface/SimpleConfigurable.h> 

#include "CLHEP/Geometry/Point3D.h"
#include "CLHEP/Units/PhysicalConstants.h"

//JetMet stuff
#include "JetMetAnalysis/MetObjects/interface/CaloEtaSliceMet.h"
#include "JetMetAnalysis/MetObjects/interface/L1CaloMet.h"
#include "JetMetAnalysis/MetObjects/interface/ConcreteMet.h"

#include "JetMetAnalysis/MetBuilders/interface/JetMetSetup.h"

#include "JetMetAnalysis/JetMetParticles/interface/JetMetParticleSetup.h"



// EtMiss class
#include <XANADOO/XANAEtMiss/interface/XANAEtMiss.h>

#include <iostream>
#include<string>

class MySetup{
public:
  MySetup();
  ~MySetup();
private:
  JetMetSetup * setup;
  double dummy;
};

MySetup::MySetup()
{
  cout<<"***********My setUp**********************";
  setup = new JetMetSetup;
}

MySetup::~MySetup()
{
  delete setup;
}

class EtMissTest : private Observer<G3EventProxy *>
{
public:
  EtMissTest();
  ~EtMissTest(){ cout<<"End of EtMissTest"<<endl;delete _mysetup;}
  void FillingTest(G3EventProxy *ev);  
private:
  void upDate(G3EventProxy *ev){if (ev!=0) FillingTest(ev);}
  void FillEtMiss(G3EventProxy *ev, EtMiss *a);
  void Print(EtMiss a);
  void PrintEvInfo(G3EventProxy *ev);
  SimpleConfigurable<string> _Scheme;
  string CheckSchemeName();
  MySetup * _mysetup;
};

EtMissTest::EtMissTest()
{    
  Observer<G3EventProxy *>::init() ;
  _Scheme=SimpleConfigurable<string>("default","EtMissScheme");
  _mysetup = new MySetup;
}

void EtMissTest::FillingTest(G3EventProxy *ev)
{
  PrintEvInfo(ev);
  EtMiss* myMet = new EtMiss(CheckSchemeName());
  FillEtMiss (ev,myMet);
  Print(*myMet);
  delete myMet;
}

void EtMissTest::PrintEvInfo(G3EventProxy *ev)
{
  static int nEvRead=0;
  cout <<"======================================================================= " << endl;
  cout << " ======> Reading " << ++nEvRead << " Event" << endl; 

  cout <<" pythia run " << ev->simSignal()->originalId().runNumber()
       <<" pythia ev  " << ev->simSignal()->originalId().eventInRun() << endl;

  cout << " oohit Run Number " << ev->simSignal()->id().runNumber()
       << " event number "     << ev->simSignal()->id().eventInRun() << endl;

  cout <<" digi run = "     << ev->pEvent()->id().runNumber()
       <<" event number = " << ev->pEvent()->id().eventInRun() << endl;

  cout <<"***************************************************************** " << endl;

}

void EtMissTest::Print(EtMiss a)
{
  cout<<"===========  scheme : "<<a.WhichScheme()<<endl;
  cout<<"*************EtMiss*************************** \n";
  cout<<" "<<a.getEMissX()<<"  "<<a.getEMissY()<<"  "<<a.getEtMiss()<<endl;
  
  cout<<"\n ********* ECAL ******************"<<endl;
  for (int i=0;i<5;i++){cout<<a.getEECal(i)<<" ";}
  cout<<endl;
  for (int i=0;i<5;i++){cout<<a.getEScalarECal(i)<<" ";}

  cout<<"\n ********* HCAL ******************"<<endl;
  for (int i=0;i<5;i++){cout<<a.getEHCal(i)<<" ";}
  cout<<endl;
  for (int i=0;i<5;i++){cout<<a.getEScalarHCal(i)<<" ";}

  cout<<"\n ********* Total ******************"<<endl;
  for (int i=0;i<5;i++){cout<<a.getE(i)<<" ";}
  cout<<endl;
  for (int i=0;i<5;i++){cout<<a.getEScalar(i)<<" ";}

  cout<<"\n ********* MC ******************"<<endl;
  for (int i=0;i<5;i++){cout<<a.getETrue(i)<<" ";}
  cout<<endl<<endl;

}

  void EtMissTest::FillEtMiss(G3EventProxy *ev, EtMiss* a)
{

  string scheme = CheckSchemeName();
    if(scheme != "default")
    {
      RecItr<ConcreteMet> gMet(ev->recEvent(),scheme.c_str());
      cout<<"AAA    ConcreteMet    AAAA"<<endl;
      while (gMet.next())
        {
          cout<<*gMet<<endl;
          a->SetEtMiss(gMet->getMissingEx(),gMet->getMissingEy());
        }
      
    }
  
  // Now it dela with ECAL and HCAL always looking at cell  
  // "vectorial"

  JetMetParticleSetup  * jps = Singleton<JetMetParticleSetup>::instance();
  double towerEtThreshold = jps->getEcalPlusHcalTowerEtCut();
  cout<<"AAA towerEtThreshold: "<<towerEtThreshold<<endl;
  double EEcal[3] = {0,0,0};
  double EHcal[3] = {0,0,0};
  double ETot[3]  = {0,0,0};
  //Scalar
  double EScEcal[5] = {0,0,0,0,0};
  double EScHcal[5] = {0,0,0,0,0};
  double EScTot[5]  = {0,0,0,0,0};

  RecItr<EcalPlusHcalTower> aTower(ev->recEvent());
  while(aTower.next())
    {
      HepPoint3D position = aTower->Position();
      double theta = position.theta();
      double phi   = position.phi();
      double EnEcal = aTower->EnergyEcalTower();
      double EnHcal = aTower->EnergyHcalTower();
      if(aTower->Energy()*sin(theta) > towerEtThreshold)
        {
          double Ex_E=EnEcal*sin(theta)*cos(phi);
          double Ey_E=EnEcal*sin(theta)*sin(phi);
          double Ez_E=EnEcal*cos(theta);
          double Ex_H=EnHcal*sin(theta)*cos(phi);
          double Ey_H=EnHcal*sin(theta)*sin(phi);
          double Ez_H=EnHcal*cos(theta);

          EEcal[0] += Ex_E;
          EEcal[1] += Ey_E;
          EEcal[2] += Ez_E;

          EHcal[0] += Ex_H;
          EHcal[1] += Ey_H;
          EHcal[2] += Ez_H;


          EScEcal[0] += fabs(Ex_E);
          EScEcal[1] += fabs(Ey_E);
          EScEcal[2] += fabs(Ez_E);
          EScEcal[3] += sqrt(Ex_E*Ex_E+Ey_E*Ey_E);
          EScEcal[4] += sqrt(Ex_E*Ex_E+Ey_E*Ey_E+Ez_E*Ez_E);

          EScHcal[0] += fabs(Ex_H);
          EScHcal[1] += fabs(Ey_H);
          EScHcal[2] += fabs(Ez_H);
          EScHcal[3] += sqrt(Ex_H*Ex_H+Ey_H*Ey_H);
          EScHcal[4] += sqrt(Ex_H*Ex_H+Ey_H*Ey_H+Ez_H*Ez_H);

        }
    }

  for (int i=0;i<3;i++)
    {
      ETot[i]=EEcal[i]+EHcal[i];
    }
  
  a->SetEECal(EEcal[0],EEcal[1],EEcal[2]);
  a->SetEHCal(EHcal[0],EHcal[1],EHcal[2]);
  
  a->SetEScalECal(EScEcal[0],EScEcal[1],EScEcal[2],EScEcal[3],EScEcal[4]);
  a->SetEScalHCal(EScHcal[0],EScHcal[1],EScHcal[2],EScHcal[3],EScHcal[4]);
  a->SetETotal(ETot[0],ETot[1],ETot[2]);
  if(scheme == "default") a->SetEtMiss(-1*ETot[0],-1*ETot[1]);
  // what EScalTotal would mean is not yet clear to me 
  a->SetEScalTotal(0.9,0.4,0.3,11,12);

}

string EtMissTest::CheckSchemeName()
{
 
  if(_Scheme.value()=="pippo") return _Scheme.value();
  if(_Scheme.value()=="pappo") return _Scheme.value();
  if(_Scheme.value()=="MetGenCalo") return _Scheme.value();
  if(_Scheme.value()=="MetCaloTower") return _Scheme.value();
  //list here all existing scheme
 
  return "default";
  
}
static PKBuilder<EtMissTest> pippo("EtMissTest") ;

---