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