TauValidation.cc

Go to the documentation of this file.

/*class TauValidation
 *  
 *  Class to fill dqm monitor elements from existing EDM file
 *
 *  $Date: 2010/09/09 11:49:20 $
 *  $Revision: 1.4 $
 */
 
#include "Validation/EventGenerator/interface/TauValidation.h"

#include "CLHEP/Units/defs.h"
#include "CLHEP/Units/PhysicalConstants.h"

#include "DataFormats/Math/interface/LorentzVector.h"

#include "DataFormats/HepMCCandidate/interface/GenParticle.h"

using namespace edm;

TauValidation::TauValidation(const edm::ParameterSet& iPSet):  
  hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")),
  tauEtCut(iPSet.getParameter<double>("tauEtCutForRtau"))
{    
  dbe = 0;
  dbe = edm::Service<DQMStore>().operator->();
}

TauValidation::~TauValidation() {}

void TauValidation::beginJob()
{
  if(dbe){
    dbe->setCurrentFolder("Generator/Tau");
    
    // Number of analyzed events
    nEvt = dbe->book1D("nEvt", "n analyzed Events", 1, 0., 1.);

    //Kinematics
    TauPt            = dbe->book1D("TauPt","Tau pT", 100 ,0,100);
    TauEta           = dbe->book1D("TauEta","Tau eta", 100 ,-2.5,2.5);
    TauProngs        = dbe->book1D("TauProngs","Tau n prongs", 7 ,0,7);
    TauDecayChannels = dbe->book1D("TauDecayChannels","Tau decay channels", 10 ,0,10);
    TauDecayChannels->setBinLabel(1+undetermined,"?");
        TauDecayChannels->setBinLabel(1+electron,"e");
        TauDecayChannels->setBinLabel(1+muon,"mu");
        TauDecayChannels->setBinLabel(1+pi,"#pi^{#pm}");
    TauDecayChannels->setBinLabel(1+pi1pi0,"#pi^{#pm}#pi^{0}");
        TauDecayChannels->setBinLabel(1+pinpi0,"#pi^{#pm}n#pi^{0}");
        TauDecayChannels->setBinLabel(1+tripi,"3#pi^{#pm}");
        TauDecayChannels->setBinLabel(1+tripinpi0,"3#pi^{#pm}n#pi^{0}");
    TauDecayChannels->setBinLabel(1+K,"K");
    TauDecayChannels->setBinLabel(1+stable,"Stable");

    TauMothers        = dbe->book1D("TauMothers","Tau mother particles", 10 ,0,10);
    TauMothers->setBinLabel(1+other,"?");
    TauMothers->setBinLabel(1+gamma,"#gamma");
    TauMothers->setBinLabel(1+Z,"Z");
    TauMothers->setBinLabel(1+W,"W");
    TauMothers->setBinLabel(1+HSM,"H_{SM}/h^{0}");
    TauMothers->setBinLabel(1+H0,"H^{0}");
    TauMothers->setBinLabel(1+A0,"A^{0}");
    TauMothers->setBinLabel(1+Hpm,"H^{#pm}");

    TauRtauW          = dbe->book1D("TauRtauW","W->Tau p(leading track)/E(visible tau)", 50 ,0,1);
    TauRtauW->setAxisTitle("rtau");
    TauRtauHpm        = dbe->book1D("TauRtauHpm","Hpm->Tau p(leading track)/E(visible tau)", 50 ,0,1);
    TauRtauHpm->setAxisTitle("rtau");
    TauSpinEffectsW   = dbe->book1D("TauSpinEffectsW","Pion energy in W rest frame", 50 ,0,1);
    TauSpinEffectsW->setAxisTitle("Energy");
    TauSpinEffectsHpm = dbe->book1D("TauSpinEffectsHpm","Pion energy in Hpm rest frame", 50 ,0,1);
    TauSpinEffectsHpm->setAxisTitle("Energy");

    TauPhotons        = dbe->book1D("TauPhoton","Photons radiating from tau", 2 ,0,2);
    TauPhotons->setBinLabel(1,"Fraction of taus radiating photons");
    TauPhotons->setBinLabel(2,"Fraction of tau pt radiated by photons");
  }

  nTaus              = 0;
  nTausWithPhotons   = 0;
  tauPtSum           = 0;
  photonFromTauPtSum = 0;

  return;
}

void TauValidation::endJob(){
  return;
}

void TauValidation::beginRun(const edm::Run& iRun,const edm::EventSetup& iSetup)
{
  iSetup.getData( fPDGTable );
  return;
}
void TauValidation::endRun(const edm::Run& iRun,const edm::EventSetup& iSetup){return;}
void TauValidation::analyze(const edm::Event& iEvent,const edm::EventSetup& iSetup)
{ 
  
  edm::Handle<HepMCProduct> evt;
  iEvent.getByLabel(hepmcCollection_, evt);

  //Get EVENT
  HepMC::GenEvent *myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));

  nEvt->Fill(0.5);

  // find taus
  for(HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); ++iter) {
    if((*iter)->status()==3) {
      if(abs((*iter)->pdg_id())==15){
        TauPt->Fill((*iter)->momentum().perp());
        TauEta->Fill((*iter)->momentum().eta());
    int mother  = tauMother(*iter);
    int decaychannel = tauDecayChannel(*iter);
    tauProngs(*iter);
    rtau(*iter,mother,decaychannel);
    spinEffects(*iter,mother,decaychannel);
    photons(*iter);
      }
    }
  }

  delete myGenEvent;
}//analyze

int TauValidation::tauMother(const HepMC::GenParticle* tau){
    int mother_pid = 0;

    if ( tau->production_vertex() ) {
        HepMC::GenVertex::particle_iterator mother;
        for (mother = tau->production_vertex()->particles_begin(HepMC::parents);
                     mother!= tau->production_vertex()->particles_end(HepMC::parents); ++mother ) {
            mother_pid = (*mother)->pdg_id();
            //std::cout << " parent " << mother_pid << std::endl;
        }
    }

    int label = other;
    if(abs(mother_pid) == 24) label = W;
        if(abs(mother_pid) == 23) label = Z;
    if(abs(mother_pid) == 22) label = gamma;
    if(abs(mother_pid) == 25) label = HSM;
    if(abs(mother_pid) == 35) label = H0;
    if(abs(mother_pid) == 36) label = A0;
    if(abs(mother_pid) == 37) label = Hpm;

    TauMothers->Fill(label);

    return mother_pid;
}

int TauValidation::tauProngs(const HepMC::GenParticle* tau){

    int nProngs = 0;
    if ( tau->end_vertex() ) {
        HepMC::GenVertex::particle_iterator des;
        for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
            des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
            int pid = (*des)->pdg_id();
            if(abs(pid) == 15) return tauProngs(*des);
            if((*des)->status() != 1) continue; // dont count unstable particles

            const HepPDT::ParticleData*  pd = fPDGTable->particle((*des)->pdg_id ());
            int charge = (int) pd->charge();
            if(charge == 0) continue;
            //std::cout << "TauValidation::tauProngs barcode=" << (*des)->barcode() << " pid=" 
                        //          << pid << " mom=" << tauMother(*des) << " status=" 
                        //          << (*des)->status() << " charge=" << charge << std::endl;
            nProngs++;
        }
    }
    TauProngs->Fill(nProngs);
    return nProngs;
}

int TauValidation::tauDecayChannel(const HepMC::GenParticle* tau){

    int channel = undetermined;
    if(tau->status() == 1) channel = stable;

    int eCount   = 0,
        muCount  = 0,
        pi0Count = 0,
            piCount  = 0;

        if ( tau->end_vertex() ) {
              HepMC::GenVertex::particle_iterator des;
              for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
                  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {

                        //if(abs(tauMother(*des)) != 15) continue;
                        int pid = (*des)->pdg_id();
                        //std::cout << " barcode=" << (*des)->barcode() << " pid="
                        //          << pid << " mom=" << tauMother(*des) << " status="
                        //          << (*des)->status() << std::endl;

                        if(abs(pid) == 15) return tauDecayChannel(*des);

            if(abs(pid) == 11)  eCount++;
            if(abs(pid) == 13)  muCount++;
            if(abs(pid) == 111) pi0Count++;
                        if(abs(pid) == 211) piCount++; 

        }
    }

    if(piCount == 1 && pi0Count == 0) channel = pi;
    if(piCount == 1 && pi0Count == 1)  channel = pi1pi0;
        if(piCount == 1 && pi0Count > 1)  channel = pinpi0;

        if(piCount == 3 && pi0Count == 0) channel = tripi;
        if(piCount == 3 && pi0Count > 0)  channel = tripinpi0;

    if(eCount == 1)                   channel = electron;
    if(muCount == 1)                  channel = muon;

    TauDecayChannels->Fill(channel);
    return channel;
}

void TauValidation::rtau(const HepMC::GenParticle* tau,int mother, int decay){

    if(decay != pi1pi0) return; // polarization only for 1-prong hadronic taus with one neutral pion to make a clean case

    if(tau->momentum().perp() < tauEtCut) return; // rtau visible only for boosted taus
    
    double rTau = 0;
    double ltrack = leadingPionMomentum(tau);
    double visibleTauE = visibleTauEnergy(tau);

    if(visibleTauE != 0) rTau = ltrack/visibleTauE;

    if(abs(mother) == 23) TauRtauW->Fill(rTau);
        if(abs(mother) == 37) TauRtauHpm->Fill(rTau); 
}

void TauValidation::spinEffects(const HepMC::GenParticle* tau,int mother, int decay){

    if(decay != pi) return; // polarization only for 1-prong hadronic taus with no neutral pions

    TLorentzVector momP4 = motherP4(tau);
    TLorentzVector pionP4 = leadingPionP4(tau);

    pionP4.Boost(-1*momP4.BoostVector());

    double energy = pionP4.E()/(momP4.M()/2);

    if(abs(mother) == 23) TauSpinEffectsW->Fill(energy);    
    if(abs(mother) == 37) TauSpinEffectsHpm->Fill(energy);
}

double TauValidation::leadingPionMomentum(const HepMC::GenParticle* tau){
    return leadingPionP4(tau).P();
}

TLorentzVector TauValidation::leadingPionP4(const HepMC::GenParticle* tau){

    TLorentzVector p4(0,0,0,0);

        if ( tau->end_vertex() ) {
              HepMC::GenVertex::particle_iterator des;
              for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
                  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
                        int pid = (*des)->pdg_id();

                        if(abs(pid) == 15) return leadingPionP4(*des);

                        if(abs(pid) != 211) continue;
 
            if((*des)->momentum().rho() > p4.P()) {
                p4 = TLorentzVector((*des)->momentum().px(),
                                                    (*des)->momentum().py(),
                                                    (*des)->momentum().pz(),
                                                    (*des)->momentum().e());
            } 
                }
        }

    return p4;
}

TLorentzVector TauValidation::motherP4(const HepMC::GenParticle* tau){

    TLorentzVector p4(0,0,0,0);

        if ( tau->production_vertex() ) {
                HepMC::GenVertex::particle_iterator mother;
                for (mother = tau->production_vertex()->particles_begin(HepMC::parents);
                     mother!= tau->production_vertex()->particles_end(HepMC::parents); ++mother ) {
                        //mother_pid = (*mother)->pdg_id();
                        //std::cout << " parent " << mother_pid << std::endl;
                        p4 = TLorentzVector((*mother)->momentum().px(),
                                            (*mother)->momentum().py(),
                                            (*mother)->momentum().pz(),
                                            (*mother)->momentum().e());
                }
        }

    return p4;
}

double TauValidation::visibleTauEnergy(const HepMC::GenParticle* tau){
    TLorentzVector p4(tau->momentum().px(),
                          tau->momentum().py(),
                          tau->momentum().pz(),
                          tau->momentum().e());

        if ( tau->end_vertex() ) {
              HepMC::GenVertex::particle_iterator des;
              for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
                  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
                        int pid = (*des)->pdg_id();

                        if(abs(pid) == 15) return visibleTauEnergy(*des);

                        if(abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
                p4 -= TLorentzVector((*des)->momentum().px(),
                                         (*des)->momentum().py(),
                                         (*des)->momentum().pz(),
                                         (*des)->momentum().e());
            }
                }
        }

    return p4.E();
}

void TauValidation::photons(const HepMC::GenParticle* tau){

        if ( tau->end_vertex() ) {
          bool photonFromTau = false;
              HepMC::GenVertex::particle_iterator des;
              for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
                  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
                        int pid = (*des)->pdg_id();
            if(pid == 22) {
                photonFromTauPtSum += (*des)->momentum().perp();
                photonFromTau = true;
            } 
              }
          nTaus++;
          if(photonFromTau) {
        tauPtSum += tau->momentum().perp();
        nTausWithPhotons++;
          }

          double nFrac = double(nTausWithPhotons)/nTaus;
          double ptFrac = 0;
          if(tauPtSum > 0) ptFrac = photonFromTauPtSum/tauPtSum;

          TauPhotons->setBinContent(1,nFrac);
          TauPhotons->setBinContent(2,ptFrac);
        }
}