Reference Manual

Go to the documentation of this file.
00001 /*class TauValidation
00002  *  
00003  *  Class to fill dqm monitor elements from existing EDM file
00004  *
00005  */
00006 #include "Validation/EventGenerator/interface/TauValidation.h"
00007 
00008 #include "CLHEP/Units/defs.h"
00009 #include "CLHEP/Units/PhysicalConstants.h"
00010 
00011 #include "DataFormats/Math/interface/LorentzVector.h"
00012 
00013 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00014 #include "Validation/EventGenerator/interface/TauDecay_CMSSW.h"
00015 #include "Validation/EventGenerator/interface/PdtPdgMini.h"
00016 
00017 using namespace edm;
00018 
00019 TauValidation::TauValidation(const edm::ParameterSet& iPSet): 
00020   _wmanager(iPSet)
00021   ,hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection"))
00022   ,tauEtCut(iPSet.getParameter<double>("tauEtCutForRtau"))
00023   ,NJAKID(22)
00024   ,zsbins(20)
00025   ,zsmin(-0.5)
00026   ,zsmax(0.5)
00027 {    
00028   dbe = 0;
00029   dbe = edm::Service<DQMStore>().operator->();
00030 }
00031 
00032 TauValidation::~TauValidation() {}
00033 
00034 void TauValidation::beginJob()
00035 {
00036   if(dbe){
00038     dbe->setCurrentFolder("Generator/Tau");
00039     
00040     // Number of analyzed events
00041     nTaus = dbe->book1D("nTaus", "n analyzed Taus", 1, 0., 1.);
00042     nPrimeTaus = dbe->book1D("nPrimeTaus", "n analyzed prime Taus", 1, 0., 1.);
00043 
00044     //Kinematics
00045     TauPt            = dbe->book1D("TauPt","Tau pT", 100 ,0,100);
00046     TauEta           = dbe->book1D("TauEta","Tau eta", 100 ,-2.5,2.5);
00047     TauPhi           = dbe->book1D("TauPhi","Tau phi", 100 ,-3.14,3.14);
00048     TauProngs        = dbe->book1D("TauProngs","Tau n prongs", 7 ,0,7);
00049     TauDecayChannels = dbe->book1D("TauDecayChannels","Tau decay channels", 13 ,0,13);
00050         TauDecayChannels->setBinLabel(1+undetermined,"?");
00051         TauDecayChannels->setBinLabel(1+electron,"e");
00052         TauDecayChannels->setBinLabel(1+muon,"mu");
00053         TauDecayChannels->setBinLabel(1+pi,"#pi^{#pm}");
00054         TauDecayChannels->setBinLabel(1+rho,"#rho^{#pm}");
00055         TauDecayChannels->setBinLabel(1+a1,"a_{1}^{#pm}");
00056         TauDecayChannels->setBinLabel(1+pi1pi0,"#pi^{#pm}#pi^{0}");
00057         TauDecayChannels->setBinLabel(1+pinpi0,"#pi^{#pm}n#pi^{0}");
00058         TauDecayChannels->setBinLabel(1+tripi,"3#pi^{#pm}");
00059         TauDecayChannels->setBinLabel(1+tripinpi0,"3#pi^{#pm}n#pi^{0}");
00060         TauDecayChannels->setBinLabel(1+K,"K");
00061         TauDecayChannels->setBinLabel(1+Kstar,"K^{*}");
00062         TauDecayChannels->setBinLabel(1+stable,"Stable");
00063 
00064     TauMothers        = dbe->book1D("TauMothers","Tau mother particles", 10 ,0,10);
00065         TauMothers->setBinLabel(1+other,"?");
00066         TauMothers->setBinLabel(1+B,"B Decays");
00067         TauMothers->setBinLabel(1+D,"D Decays");
00068         TauMothers->setBinLabel(1+gamma,"#gamma");
00069         TauMothers->setBinLabel(1+Z,"Z");
00070         TauMothers->setBinLabel(1+W,"W");
00071         TauMothers->setBinLabel(1+HSM,"H_{SM}/h^{0}");
00072         TauMothers->setBinLabel(1+H0,"H^{0}");
00073         TauMothers->setBinLabel(1+A0,"A^{0}");
00074         TauMothers->setBinLabel(1+Hpm,"H^{#pm}");
00075 
00076     TauRtauW          = dbe->book1D("TauRtauW","W->Tau p(leading track)/E(visible tau)", 50 ,0,1);     TauRtauW->setAxisTitle("rtau");
00077     TauRtauHpm        = dbe->book1D("TauRtauHpm","Hpm->Tau p(leading track)/E(visible tau)", 50 ,0,1); TauRtauHpm->setAxisTitle("rtau");
00078 
00079     TauSpinEffectsW_X   = dbe->book1D("TauSpinEffectsWX","Pion energy in W rest frame", 50 ,0,1);     TauSpinEffectsW_X->setAxisTitle("X");
00080     TauSpinEffectsHpm_X = dbe->book1D("TauSpinEffectsHpmX","Pion energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_X->setAxisTitle("X");
00081 
00082     TauSpinEffectsW_eX   = dbe->book1D("TauSpinEffectsWeX","e energy in W rest frame", 50 ,0,1);     TauSpinEffectsW_eX->setAxisTitle("X");
00083     TauSpinEffectsHpm_eX = dbe->book1D("TauSpinEffectsHpmeX","e energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_eX->setAxisTitle("X");
00084 
00085     TauSpinEffectsW_muX   = dbe->book1D("TauSpinEffectsWmuX","mu energy in W rest frame", 50 ,0,1);     TauSpinEffectsW_muX->setAxisTitle("X");
00086     TauSpinEffectsHpm_muX = dbe->book1D("TauSpinEffectsHpmmuX","mu energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_muX->setAxisTitle("X");
00087 
00088     TauSpinEffectsW_UpsilonRho   = dbe->book1D("TauSpinEffectsWUpsilonRho","#Upsilon for #rho", 50 ,-1,1);     TauSpinEffectsW_UpsilonRho->setAxisTitle("#Upsilon");
00089     TauSpinEffectsHpm_UpsilonRho = dbe->book1D("TauSpinEffectsHpmUpsilonRho","#Upsilon for #rho", 50 ,-1,1);   TauSpinEffectsHpm_UpsilonRho->setAxisTitle("#Upsilon");
00090 
00091     TauSpinEffectsW_UpsilonA1   = dbe->book1D("TauSpinEffectsWUpsilonA1","#Upsilon for a1", 50 ,-1,1);       TauSpinEffectsW_UpsilonA1->setAxisTitle("#Upsilon");
00092     TauSpinEffectsHpm_UpsilonA1 = dbe->book1D("TauSpinEffectsHpmUpsilonA1","#Upsilon for a1", 50 ,-1,1);     TauSpinEffectsHpm_UpsilonA1->setAxisTitle("#Upsilon");
00093 
00094     TauSpinEffectsZ_MVis   = dbe->book1D("TauSpinEffectsZMVis","Mass of pi+ pi-", 25 ,0,1.1);       TauSpinEffectsZ_MVis->setAxisTitle("M_{#pi^{+}#pi^{-}}");
00095     TauSpinEffectsH_MVis   = dbe->book1D("TauSpinEffectsHMVis","Mass of pi+ pi-", 25 ,0,1.1);       TauSpinEffectsZ_MVis->setAxisTitle("M_{#pi^{+}#pi^{-}}");
00096 
00097     TauSpinEffectsZ_Zs   = dbe->book1D("TauSpinEffectsZZs","Z_{s}", zsbins ,zsmin,zsmax);        TauSpinEffectsZ_Zs->setAxisTitle("Z_{s}");
00098     TauSpinEffectsH_Zs   = dbe->book1D("TauSpinEffectsHZs","Z_{s}", zsbins ,zsmin,zsmax);        TauSpinEffectsZ_Zs->setAxisTitle("Z_{s}");
00099 
00100     TauSpinEffectsZ_Xf   = dbe->book1D("TauSpinEffectsZXf","X of forward emitted #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_Xf->setAxisTitle("X_{f}");
00101     TauSpinEffectsH_Xf   = dbe->book1D("TauSpinEffectsHXf","X of forward emitted #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_Xf->setAxisTitle("X_{f}");
00102 
00103     TauSpinEffectsZ_Xb   = dbe->book1D("TauSpinEffectsZXb","X of backward emitted #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_Xb->setAxisTitle("X_{b}");
00104     TauSpinEffectsH_Xb   = dbe->book1D("TauSpinEffectsHXb","X of backward emitted #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_Xb->setAxisTitle("X_{b}");
00105 
00106     TauSpinEffectsZ_eX   = dbe->book1D("TauSpinEffectsZeX","e energy in Z rest frame", 50 ,0,1);     TauSpinEffectsZ_eX->setAxisTitle("X");
00107     TauSpinEffectsH_eX = dbe->book1D("TauSpinEffectsHeX","e energy in H rest frame", 50 ,0,1); TauSpinEffectsH_eX->setAxisTitle("X");
00108 
00109     TauSpinEffectsZ_muX   = dbe->book1D("TauSpinEffectsZmuX","mu energy in z rest frame", 50 ,0,1);     TauSpinEffectsZ_muX->setAxisTitle("X");
00110     TauSpinEffectsH_muX = dbe->book1D("TauSpinEffectsHmuX","mu energy in H rest frame", 50 ,0,1); TauSpinEffectsH_muX->setAxisTitle("X");
00111 
00112 
00113     TauFSRPhotonsN=dbe->book1D("TauFSRPhotonsN","FSR Photons radiating from/with tau (Gauge Boson)", 5 ,-0.5,4.5);
00114     TauFSRPhotonsN->setAxisTitle("N FSR Photons radiating from/with tau");
00115     TauFSRPhotonsPt=dbe->book1D("TauFSRPhotonsPt","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
00116     TauFSRPhotonsPt->setAxisTitle("P_{t} of FSR Photons radiating from/with tau [per tau]");
00117     TauFSRPhotonsPtSum=dbe->book1D("TauFSRPhotonsPtSum","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
00118     TauFSRPhotonsPtSum->setAxisTitle("P_{t} of FSR Photons radiating from/with tau [per tau]");
00119 
00120     TauBremPhotonsN=dbe->book1D("TauBremPhotonsN","Brem. Photons radiating in tau decay", 5 ,-0.5,4.5);
00121     TauBremPhotonsN->setAxisTitle("N FSR Photons radiating from/with tau");
00122     TauBremPhotonsPt=dbe->book1D("TauBremPhotonsPt","Sum Brem Pt ", 100 ,0,100);
00123     TauFSRPhotonsPt->setAxisTitle("P_{t} of Brem. Photons radiating in tau decay");    
00124     TauBremPhotonsPtSum =dbe->book1D("TauBremPhotonsPtSum","Sum of Brem Pt ", 100 ,0,100);
00125     TauFSRPhotonsPtSum->setAxisTitle("Sum P_{t} of Brem. Photons radiating in tau decay");
00126 
00127     JAKID =dbe->book1D("JAKID","JAK ID",NJAKID+1,-0.5,NJAKID+0.5);
00128     for(unsigned int i=0; i<NJAKID+1;i++){
00129       JAKInvMass.push_back(std::vector<MonitorElement *>());
00130       TString tmp="JAKID";
00131       tmp+=i;
00132       JAKInvMass.at(i).push_back(dbe->book1D("M"+tmp,"M_{"+tmp+"} (GeV)", 80 ,0,2.0));
00133       if(i==TauDecay::JAK_A1_3PI ||
00134          i==TauDecay::JAK_KPIK ||
00135          i==TauDecay::JAK_KPIPI ){
00136         JAKInvMass.at(i).push_back(dbe->book1D("M13"+tmp,"M_{13,"+tmp+"} (GeV)", 80 ,0,2.0));
00137         JAKInvMass.at(i).push_back(dbe->book1D("M23"+tmp,"M_{23,"+tmp+"} (GeV)", 80 ,0,2.0));
00138         JAKInvMass.at(i).push_back(dbe->book1D("M12"+tmp,"M_{12,"+tmp+"} (GeV)", 80 ,0,2.0));
00139       }
00140     }
00141   }
00142   
00143   return;
00144 }
00145 
00146 void TauValidation::endJob(){
00147   return;
00148 }
00149 
00150 void TauValidation::beginRun(const edm::Run& iRun,const edm::EventSetup& iSetup)
00151 {
00153   iSetup.getData( fPDGTable );
00154   return;
00155 }
00156 void TauValidation::endRun(const edm::Run& iRun,const edm::EventSetup& iSetup){return;}
00157 void TauValidation::analyze(const edm::Event& iEvent,const edm::EventSetup& iSetup)
00158 { 
00160   edm::Handle<HepMCProduct> evt;
00161   iEvent.getByLabel(hepmcCollection_, evt);
00162 
00163   //Get EVENT
00164   HepMC::GenEvent *myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
00165 
00166   double weight = _wmanager.weight(iEvent);
00167 
00169   /*
00170   edm::Handle<double> WT;
00171   iEvent.getByLabel(edm::InputTag("TauSpinnerGen","TauSpinnerWT"),WT);
00172   weight = 1.0;
00173   if(*(WT.product())>1e-3 && *(WT.product())<=10.0) weight=(*(WT.product()));
00174   else {weight=1.0;}
00175   */
00177 
00178   // find taus
00179   for(HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); iter++) {
00180     if(abs((*iter)->pdg_id())==23){
00181       spinEffectsZ(*iter,weight);
00182     }
00183     if(abs((*iter)->pdg_id())==15){
00184       if(isLastTauinChain(*iter)){
00185         nTaus->Fill(0.5,weight);
00186         int mother  = tauMother(*iter,weight);
00187         int decaychannel = tauDecayChannel(*iter,weight);
00188         tauProngs(*iter, weight);
00189         if(mother>-1){ // exclude B, D and other non-signal decay modes
00190           nPrimeTaus->Fill(0.5,weight);
00191           TauPt->Fill((*iter)->momentum().perp(),weight);
00192           TauEta->Fill((*iter)->momentum().eta(),weight);
00193           TauPhi->Fill((*iter)->momentum().phi(),weight);
00194           rtau(*iter,mother,decaychannel,weight);
00195           photons(*iter,weight);
00196         }
00198         //Adding JAKID and Mass information
00199         //
00200         TauDecay_CMSSW TD;
00201         unsigned int jak_id, TauBitMask;
00202         if(TD.AnalyzeTau((*iter),jak_id,TauBitMask,false,false)){
00203           JAKID->Fill(jak_id,weight);
00204           if(jak_id<=NJAKID){
00205             int tcharge=(*iter)->pdg_id()/abs((*iter)->pdg_id());
00206             std::vector<HepMC::GenParticle*> part=TD.Get_TauDecayProducts();
00207             spinEffects(*iter,mother,jak_id,part,weight);
00208             TLorentzVector LVQ(0,0,0,0);
00209             TLorentzVector LVS12(0,0,0,0);
00210             TLorentzVector LVS13(0,0,0,0);
00211             TLorentzVector LVS23(0,0,0,0);
00212             bool haspart1=false;
00213             for(unsigned int i=0;i<part.size();i++){
00214               if(TD.isTauFinalStateParticle(part.at(i)->pdg_id()) &&
00215                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_e &&
00216                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_mu &&
00217                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_tau ){
00218                 TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00219                 LVQ+=LV;
00220                 if(jak_id==TauDecay::JAK_A1_3PI ||
00221                    jak_id==TauDecay::JAK_KPIK ||
00222                    jak_id==TauDecay::JAK_KPIPI
00223                    ){
00224                   if((tcharge==part.at(i)->pdg_id()/abs(part.at(i)->pdg_id()) && TD.nProng(TauBitMask)==3) || (jak_id==TauDecay::JAK_A1_3PI && TD.nProng(TauBitMask)==1 && abs(part.at(i)->pdg_id())==PdtPdgMini::pi_plus) ){
00225                     LVS13+=LV;
00226                     LVS23+=LV;
00227                   }
00228                   else{
00229                     LVS12+=LV;
00230                     if(!haspart1 && ((jak_id==TauDecay::JAK_A1_3PI)  || (jak_id!=TauDecay::JAK_A1_3PI && abs(part.at(i)->pdg_id())==PdtPdgMini::K_plus) )){
00231                       LVS13+=LV;
00232                       haspart1=true;
00233                     }
00234                     else{
00235                       LVS23+=LV;
00236                     }
00237                   }
00238                 }
00239               }
00240             }
00241             part.clear();
00242             JAKInvMass.at(jak_id).at(0)->Fill(LVQ.M(),weight);
00243             if(jak_id==TauDecay::JAK_A1_3PI ||
00244                jak_id==TauDecay::JAK_KPIK ||
00245                jak_id==TauDecay::JAK_KPIPI
00246                ){
00247               JAKInvMass.at(jak_id).at(1)->Fill(LVS13.M(),weight);
00248               JAKInvMass.at(jak_id).at(2)->Fill(LVS23.M(),weight);
00249               JAKInvMass.at(jak_id).at(3)->Fill(LVS12.M(),weight);
00250             }
00251           }
00252         }
00253         else{
00254           JAKID->Fill(jak_id,weight);  
00255         }
00256       }
00257     }
00258   }
00259   delete myGenEvent;
00260 }//analyze
00261 
00262 const HepMC::GenParticle* TauValidation::GetMother(const HepMC::GenParticle* tau){
00263   if ( tau->production_vertex() ) {
00264     HepMC::GenVertex::particle_iterator mother;
00265     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00266       if((*mother)->pdg_id() == tau->pdg_id()) return GetMother(*mother);
00267       return (*mother);
00268     }
00269   }
00270   return tau;
00271 }
00272 
00273 
00274 const std::vector<HepMC::GenParticle*> TauValidation::GetMothers(const HepMC::GenParticle* boson){
00275   std::vector<HepMC::GenParticle*> mothers;
00276   if ( boson->production_vertex() ) {
00277     HepMC::GenVertex::particle_iterator mother;
00278     for (mother = boson->production_vertex()->particles_begin(HepMC::parents); mother!= boson->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00279       if((*mother)->pdg_id() == boson->pdg_id()) return GetMothers(*mother);
00280       mothers.push_back(*mother);
00281     }
00282   }
00283   return mothers;
00284 }
00285 
00286 
00287 int TauValidation::findMother(const HepMC::GenParticle* tau){
00288   int mother_pid = 0;
00289   if ( tau->production_vertex() ) {
00290     HepMC::GenVertex::particle_iterator mother;
00291     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00292       mother_pid = (*mother)->pdg_id();
00293       if(mother_pid == tau->pdg_id()) return findMother(*mother); //mother_pid = -1; Make recursive to look for last tau in chain
00294     }
00295   }
00296   return mother_pid;
00297 }
00298 
00299 
00300 bool TauValidation::isLastTauinChain(const HepMC::GenParticle* tau){
00301   if ( tau->end_vertex() ) {
00302     HepMC::GenVertex::particle_iterator dau;
00303     for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
00304       int dau_pid = (*dau)->pdg_id();
00305       if(dau_pid == tau->pdg_id()) return false;
00306     }
00307   }
00308   return true;
00309 }
00310 
00311 
00312 void TauValidation::findTauList(const HepMC::GenParticle* tau,std::vector<const HepMC::GenParticle*> &TauList){
00313   TauList.insert(TauList.begin(),tau);
00314   if ( tau->production_vertex() ) {
00315     HepMC::GenVertex::particle_iterator mother;
00316     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents);mother++) {
00317       if((*mother)->pdg_id() == tau->pdg_id()){
00318         findTauList(*mother,TauList);
00319       }
00320     }
00321   }
00322 }
00323 
00324 void TauValidation::findFSRandBrem(const HepMC::GenParticle* p, bool doBrem, std::vector<const HepMC::GenParticle*> &ListofFSR,
00325                                   std::vector<const HepMC::GenParticle*> &ListofBrem){
00326   // note this code split the FSR and Brem based one if the tau decays into a tau+photon or not with the Fortran Tauola Interface, this is not 100% correct because photos puts the tau with the regular tau decay products. 
00327   if(abs(p->pdg_id())==15){
00328     if(isLastTauinChain(p)){ doBrem=true;}
00329     else{ doBrem=false;}
00330   }
00331     int photo_ID=22;
00332   if ( p->end_vertex() ) {
00333     HepMC::GenVertex::particle_iterator dau;
00334     for (dau = p->end_vertex()->particles_begin(HepMC::children); dau!= p->end_vertex()->particles_end(HepMC::children); dau++ ) {
00335       //if(doBrem) std::cout << "true " << (*dau)->pdg_id() << " "  << std::endl;
00336       //else std::cout << "false " << (*dau)->pdg_id() << " " << std::endl;
00337       if(abs((*dau)->pdg_id()) == abs(photo_ID) && !doBrem){ListofFSR.push_back(*dau);}
00338       if(abs((*dau)->pdg_id()) == abs(photo_ID) && doBrem){ListofBrem.push_back(*dau);}
00339       if((*dau)->end_vertex() && (*dau)->end_vertex()->particles_out_size()>0 && abs((*dau)->pdg_id()) != 111 && abs((*dau)->pdg_id()) != 221/* remove pi0 and eta decays*/){
00340         findFSRandBrem(*dau,doBrem,ListofFSR,ListofBrem);
00341       }
00342     }
00343   }
00344 }
00345 
00346 
00347 
00348 void TauValidation::FindPhotosFSR(const HepMC::GenParticle* p,std::vector<const HepMC::GenParticle*> &ListofFSR,double &BosonScale){
00349   BosonScale=0.0;
00350   const HepMC::GenParticle* m=GetMother(p);
00351   double mother_pid=m->pdg_id();
00352   if(m->end_vertex() && mother_pid!=p->pdg_id()){
00353     HepMC::GenVertex::particle_iterator dau;
00354     for (dau = m->end_vertex()->particles_begin(HepMC::children); dau!= m->end_vertex()->particles_end(HepMC::children); dau++ ) {
00355       int dau_pid = (*dau)->pdg_id();
00356       if(fabs(dau_pid) == 22) {
00357         ListofFSR.push_back(*dau);
00358       }
00359     }
00360   }
00361   if(abs(mother_pid) == 24) BosonScale=1.0; // W
00362   if(abs(mother_pid) == 23) BosonScale=2.0; // Z;
00363   if(abs(mother_pid) == 22) BosonScale=2.0; // gamma;
00364   if(abs(mother_pid) == 25) BosonScale=2.0; // HSM;
00365   if(abs(mother_pid) == 35) BosonScale=2.0; // H0;
00366   if(abs(mother_pid) == 36) BosonScale=2.0; // A0;
00367   if(abs(mother_pid) == 37) BosonScale=1.0; //Hpm;
00368 }
00369 
00370 
00371 int TauValidation::tauMother(const HepMC::GenParticle* tau, double weight){
00372 
00373   if(abs(tau->pdg_id()) != 15 ) return -3;
00374   
00375   int mother_pid = findMother(tau);
00376   if(mother_pid == -2) return -2;
00377   
00378   int label = other;
00379   if(abs(mother_pid) == 24) label = W;
00380   if(abs(mother_pid) == 23) label = Z;
00381   if(abs(mother_pid) == 22) label = gamma;
00382   if(abs(mother_pid) == 25) label = HSM;
00383   if(abs(mother_pid) == 35) label = H0;
00384   if(abs(mother_pid) == 36) label = A0;
00385   if(abs(mother_pid) == 37) label = Hpm;
00386   
00387   int mother_shortpid=(abs(mother_pid)%10000);
00388   if(mother_shortpid>500 && mother_shortpid<600 )label = B;
00389   if(mother_shortpid>400 && mother_shortpid<500)label = D;
00390   TauMothers->Fill(label,weight);
00391   if(label==B || label == D || label == other) return -1;
00392   
00393   return mother_pid;
00394 }
00395 
00396 int TauValidation::tauProngs(const HepMC::GenParticle* tau, double weight){
00397   int nProngs = 0;
00398   if ( tau->end_vertex() ) {
00399     HepMC::GenVertex::particle_iterator des;
00400     for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00401         des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00402       int pid = (*des)->pdg_id();
00403       if(abs(pid) == 15) return tauProngs(*des, weight);
00404       if((*des)->status() != 1) continue; // dont count unstable particles
00405       
00406       const HepPDT::ParticleData*  pd = fPDGTable->particle((*des)->pdg_id ());
00407       int charge = (int) pd->charge();
00408       if(charge == 0) continue;
00409       nProngs++;
00410     }
00411   }
00412   TauProngs->Fill(nProngs,weight);
00413   return nProngs;
00414 }
00415 
00416 int TauValidation::tauDecayChannel(const HepMC::GenParticle* tau, double weight){
00417 
00418   int channel = undetermined;
00419   if(tau->status() == 1) channel = stable;
00420   int allCount   = 0,
00421     eCount     = 0,
00422     muCount    = 0,
00423     pi0Count   = 0,
00424     piCount    = 0,
00425     rhoCount   = 0,
00426     a1Count    = 0,
00427     KCount     = 0,
00428     KstarCount = 0;
00429   
00430   if ( tau->end_vertex() ) {
00431     HepMC::GenVertex::particle_iterator des;
00432     for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00433         des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00434       int pid = (*des)->pdg_id();
00435       if(abs(pid) == 15) return tauDecayChannel(*des,weight);
00436       
00437       allCount++;
00438       if(abs(pid) == 11)    eCount++;
00439       if(abs(pid) == 13)    muCount++;
00440       if(abs(pid) == 111)   pi0Count++;
00441       if(abs(pid) == 211)   piCount++;
00442       if(abs(pid) == 213)   rhoCount++;
00443       if(abs(pid) == 20213) a1Count++;
00444       if(abs(pid) == 321)   KCount++;
00445       if(abs(pid) == 323)   KstarCount++;
00446       
00447     }
00448   }
00449   // resonances  
00450   if(KCount >= 1)     channel = K;
00451   if(KstarCount >= 1) channel = Kstar;
00452   if(a1Count >= 1)    channel = a1;
00453   if(rhoCount >= 1)   channel = rho;
00454   if(channel!=undetermined && weight!=0.0) TauDecayChannels->Fill(channel,weight);
00455   
00456   // final state products
00457   if(piCount == 1 && pi0Count == 0) channel = pi;
00458   if(piCount == 1 && pi0Count == 1) channel = pi1pi0;
00459   if(piCount == 1 && pi0Count > 1)  channel = pinpi0;
00460   if(piCount == 3 && pi0Count == 0) channel = tripi;
00461   if(piCount == 3 && pi0Count > 0)  channel = tripinpi0;
00462   if(eCount == 1)                   channel = electron;
00463   if(muCount == 1)                  channel = muon;
00464   if(weight!=0.0) TauDecayChannels->Fill(channel,weight);
00465   return channel;
00466 }
00467 
00468 
00469 void TauValidation::rtau(const HepMC::GenParticle* tau,int mother, int decay, double weight){
00470 
00471         if(decay != pi1pi0) return; // polarization only for 1-prong hadronic taus with one neutral pion to make a clean case
00472 
00473         if(tau->momentum().perp() < tauEtCut) return; // rtau visible only for boosted taus
00474         
00475         double rTau = 0;
00476         double ltrack = leadingPionMomentum(tau, weight);
00477         double visibleTauE = visibleTauEnergy(tau);
00478 
00479         if(visibleTauE != 0) rTau = ltrack/visibleTauE;
00480 
00481         if(abs(mother) == 24) TauRtauW->Fill(rTau,weight);
00482         if(abs(mother) == 37) TauRtauHpm->Fill(rTau,weight); 
00483 }
00484 
00485 void TauValidation::spinEffects(const HepMC::GenParticle* tau,int mother, int decay, std::vector<HepMC::GenParticle*> &part,double weight){
00486   if(decay == TauDecay::JAK_PION || decay == TauDecay::JAK_MUON || decay == TauDecay::JAK_ELECTRON){  // polarization only for 1-prong hadronic taus with no neutral pions
00487     TLorentzVector momP4 = motherP4(tau);
00488     TLorentzVector pionP4 = leadingPionP4(tau);
00489     pionP4.Boost(-1*momP4.BoostVector());
00490     double energy = pionP4.E()/(momP4.M()/2);
00491     if(decay == TauDecay::JAK_PION){
00492       if(abs(mother) == 24) TauSpinEffectsW_X->Fill(energy,weight);     
00493       if(abs(mother) == 37) TauSpinEffectsHpm_X->Fill(energy,weight);
00494     }
00495     if(decay == TauDecay::JAK_MUON){
00496       if(abs(mother) == 24) TauSpinEffectsW_muX->Fill(energy,weight);
00497       if(abs(mother) == 37) TauSpinEffectsHpm_muX->Fill(energy,weight);
00498     }
00499     if(decay == TauDecay::JAK_ELECTRON){
00500       if(abs(mother) == 24) TauSpinEffectsW_eX->Fill(energy,weight);
00501       if(abs(mother) == 37) TauSpinEffectsHpm_eX->Fill(energy,weight);
00502     }
00503 
00504   }
00505   else if(decay==TauDecay::JAK_RHO_PIPI0){
00506     TLorentzVector rho(0,0,0,0),pi(0,0,0,0);
00507     for(unsigned int i=0;i<part.size();i++){
00508       TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00509       if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi_plus){pi+=LV; rho+=LV;}
00510       if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi0){rho+=LV;}
00511     }
00512     if(abs(mother) == 24) TauSpinEffectsW_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
00513     if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
00514   }
00515   else if(decay==TauDecay::JAK_A1_3PI){ // only for pi2pi0 for now
00516     TLorentzVector a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0); 
00517     int nplus(0),nminus(0);
00518     for(unsigned int i=0;i<part.size();i++){
00519       TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00520       if(part.at(i)->pdg_id()==PdtPdgMini::pi_plus){ pi_p+=LV; a1+=LV; nplus++;}
00521       if(part.at(i)->pdg_id()==PdtPdgMini::pi_minus){pi_m+=LV; a1+=LV; nminus++;}
00522     }
00523     double gamma=0;
00524     if(nplus+nminus==3 && nplus==1)  gamma=2*pi_p.Pt()/a1.Pt()-1;
00525     if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.Pt()/a1.Pt()-1;
00526     else{
00527       pi_p+=pi_m; gamma=2*pi_p.Pt()/a1.Pt()-1;
00528     }
00529     if(abs(mother) == 24) TauSpinEffectsW_UpsilonA1->Fill(gamma,weight);
00530     if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonA1->Fill(gamma,weight);
00531   }
00532 }
00533 
00534 void TauValidation::spinEffectsZ(const HepMC::GenParticle* boson, double weight){
00535 
00536   TLorentzVector tautau(0,0,0,0);
00537   TLorentzVector pipi(0,0,0,0);
00538   TLorentzVector taum(0,0,0,0);
00539   int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
00540   double x1(0),x2(0); 
00541   TLorentzVector Zboson(boson->momentum().px(),boson->momentum().py(),boson->momentum().pz(),boson->momentum().e());
00542   if ( boson->end_vertex() ) {
00543     HepMC::GenVertex::particle_iterator des;
00544     for(des = boson->end_vertex()->particles_begin(HepMC::children);    des!= boson->end_vertex()->particles_end(HepMC::children);++des ) {
00545       int pid = (*des)->pdg_id();
00546       if(abs(findMother(*des)) != 15 &&  abs(pid) == 15 && (tauDecayChannel(*des) == pi || tauDecayChannel(*des) == muon || tauDecayChannel(*des) == electron )){
00547         if(tauDecayChannel(*des) == pi)nSinglePionDecays++;
00548         if(tauDecayChannel(*des) == muon)nSingleMuonDecays++;
00549         if(tauDecayChannel(*des) == electron)nSingleElectronDecays++;
00550         TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
00551         tautau += LVtau;
00552         TLorentzVector LVpi=leadingPionP4(*des);
00553         pipi+=LVpi;
00554         const HepPDT::ParticleData*  pd = fPDGTable->particle((*des)->pdg_id ());
00555         int charge = (int) pd->charge();
00556         LVtau.Boost(-1*Zboson.BoostVector());
00557         LVpi.Boost(-1*Zboson.BoostVector());
00558         if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
00559         else{ x2=LVpi.P()/LVtau.E();}
00560      }
00561     }
00562   }
00563   if(nSingleMuonDecays==2){
00564     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_muX->Fill(x1,weight);
00565     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_muX->Fill(x1,weight);
00566   }
00567   if(nSingleElectronDecays==2){
00568     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_eX->Fill(x1,weight);
00569     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_eX->Fill(x1,weight);
00570   }
00571   if(nSinglePionDecays == 2 && tautau.M()!= 0) {
00572     for(int i=0;i<zsbins;i++){
00573       double zslow=((double)i)*(zsmax-zsmin)/((double)zsbins)+zsmin; 
00574       double zsup=((double)i+1)*(zsmax-zsmin)/((double)zsbins)+zsmin;
00575       double aup=Zstoa(zsup), alow=Zstoa(zslow);
00576       if(x2-x1>alow && x2-x1<aup){
00577         double zs=(zsup+zslow)/2;
00578         if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_Zs->Fill(zs,weight);
00579         if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Zs->Fill(zs,weight);
00580         break;
00581       }
00582     }
00583     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_MVis->Fill(pipi.M()/tautau.M(),weight);
00584     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_MVis->Fill(pipi.M()/tautau.M(),weight);
00585  
00586     if(x1!=0){
00587       const std::vector<HepMC::GenParticle*> m=GetMothers(boson);
00588       int q(0),qbar(0);
00589       TLorentzVector Z(0,0,0,0);
00590       for(unsigned int i=0;i<m.size();i++){
00591         if(m.at(i)->pdg_id()==PdtPdgMini::d      || m.at(i)->pdg_id()==PdtPdgMini::u      ){q++;}
00592         if(m.at(i)->pdg_id()==PdtPdgMini::anti_d || m.at(i)->pdg_id()==PdtPdgMini::anti_u ){qbar++;}
00593       }
00594       if(q==1 && qbar==1){// assume q has largest E (valence vs see quarks) 
00595         if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
00596           if(abs(boson->pdg_id())==PdtPdgMini::Z0)      TauSpinEffectsZ_Xf->Fill(x1,weight);
00597           if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xf->Fill(x1,weight);
00598         }
00599         else{
00600           if(abs(boson->pdg_id())==PdtPdgMini::Z0)      TauSpinEffectsZ_Xb->Fill(x1,weight);
00601           if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xb->Fill(x1,weight);
00602         }
00603       }
00604     }
00605   }
00606 }
00607 
00608 double TauValidation::Zstoa(double zs){
00609   double a=1-sqrt(fabs(1.0-2*fabs(zs)));
00610   if(zs<0){
00611     a*=-1.0;
00612   }
00613   return a;
00614 }
00615 
00616 
00617 double TauValidation::leadingPionMomentum(const HepMC::GenParticle* tau, double weight){
00618         return leadingPionP4(tau).P();
00619 }
00620 
00621 TLorentzVector TauValidation::leadingPionP4(const HepMC::GenParticle* tau){
00622 
00623         TLorentzVector p4(0,0,0,0);
00624 
00625         if ( tau->end_vertex() ) {
00626               HepMC::GenVertex::particle_iterator des;
00627               for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00628                   des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00629                         int pid = (*des)->pdg_id();
00630 
00631                         if(abs(pid) == 15) return leadingPionP4(*des);
00632 
00633                         if(!(abs(pid)==211 || abs(pid)==13 || abs(pid)==11)) continue;
00634  
00635                         if((*des)->momentum().rho() > p4.P()) {
00636                                 p4 = TLorentzVector((*des)->momentum().px(),
00637                                                     (*des)->momentum().py(),
00638                                                     (*des)->momentum().pz(),
00639                                                     (*des)->momentum().e());
00640                         } 
00641                 }
00642         }
00643 
00644         return p4;
00645 }
00646 
00647 TLorentzVector TauValidation::motherP4(const HepMC::GenParticle* tau){
00648   const HepMC::GenParticle* m=GetMother(tau);
00649   return TLorentzVector(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
00650 }
00651 
00652 double TauValidation::visibleTauEnergy(const HepMC::GenParticle* tau){
00653         TLorentzVector p4(tau->momentum().px(),
00654                           tau->momentum().py(),
00655                           tau->momentum().pz(),
00656                           tau->momentum().e());
00657 
00658         if ( tau->end_vertex() ) {
00659               HepMC::GenVertex::particle_iterator des;
00660               for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00661                   des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00662                         int pid = (*des)->pdg_id();
00663 
00664                         if(abs(pid) == 15) return visibleTauEnergy(*des);
00665 
00666                         if(abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
00667                                 p4 -= TLorentzVector((*des)->momentum().px(),
00668                                                      (*des)->momentum().py(),
00669                                                      (*des)->momentum().pz(),
00670                                                      (*des)->momentum().e());
00671                         }
00672                 }
00673         }
00674 
00675         return p4.E();
00676 }
00677 
00678 void TauValidation::photons(const HepMC::GenParticle* tau, double weight){
00679   // Find First tau in chain
00680   std::vector<const HepMC::GenParticle*> TauList;
00681   findTauList(tau,TauList);
00682 
00683   // Get List of Gauge Boson to tau(s) FSR and Brem
00684   bool passedW=false;
00685   std::vector<const HepMC::GenParticle*> ListofFSR;  ListofFSR.clear();
00686   std::vector<const HepMC::GenParticle*> ListofBrem; ListofBrem.clear();
00687   std::vector<const HepMC::GenParticle*> FSR_photos; FSR_photos.clear();
00688   double BosonScale(1);
00689   if(TauList.size()>0){
00690     TauValidation::findFSRandBrem(TauList.at(0),passedW,ListofFSR,ListofBrem);
00691     TauValidation::FindPhotosFSR(TauList.at(0),FSR_photos,BosonScale);
00692 
00693     // Add the Tau Brem. information
00694     TauBremPhotonsN->Fill(ListofBrem.size(),weight);
00695     double photonPtSum=0;
00696     for(unsigned int i=0;i<ListofBrem.size();i++){
00697       photonPtSum+=ListofBrem.at(i)->momentum().perp();
00698       TauBremPhotonsPt->Fill(ListofBrem.at(i)->momentum().perp(),weight);
00699     }
00700     TauBremPhotonsPtSum->Fill(photonPtSum,weight);
00701         
00702     // Now add the Gauge Boson FSR information
00703     if(BosonScale!=0){
00704       TauFSRPhotonsN->Fill(ListofFSR.size(),weight);
00705       photonPtSum=0;
00706       for(unsigned int i=0;i<ListofFSR.size();i++){
00707         photonPtSum+=ListofFSR.at(i)->momentum().perp();
00708         TauFSRPhotonsPt->Fill(ListofFSR.at(i)->momentum().perp(),weight);
00709       }
00710       double FSR_photosSum(0);
00711       for(unsigned int i=0;i<FSR_photos.size();i++){
00712         FSR_photosSum+=FSR_photos.at(i)->momentum().perp();
00713         TauFSRPhotonsPt->Fill(FSR_photos.at(i)->momentum().perp()/BosonScale,weight*BosonScale);
00714       }
00715       TauFSRPhotonsPtSum->Fill(photonPtSum+FSR_photosSum/BosonScale,weight);
00716     }
00717   }
00718 }
00719
/data/git/CMSSW_5_3_11_patch5/src/Validation/EventGenerator/plugins/TauValidation.cc
