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_X= dbe->book1D("TauSpinEffectsZX","X of #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_X->setAxisTitle("X");
00101     TauSpinEffectsH_X= dbe->book1D("TauSpinEffectsH_X","X of #tau^{-}", 25 ,0,1.0);           TauSpinEffectsH_X->setAxisTitle("X");
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 { 
00163   edm::Handle<HepMCProduct> evt;
00164   iEvent.getByLabel(hepmcCollection_, evt);
00165 
00166   //Get EVENT
00167   HepMC::GenEvent *myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
00168 
00169   double weight = _wmanager.weight(iEvent);
00170 
00172   /*
00173   edm::Handle<double> WT;
00174   iEvent.getByLabel(edm::InputTag("TauSpinnerGen","TauSpinnerWT"),WT);
00175   weight = 1.0;
00176   if(*(WT.product())>1e-3 && *(WT.product())<=10.0) weight=(*(WT.product()));
00177   else {weight=1.0;}
00178   */
00180 
00181   // find taus
00182   for(HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); iter++) {
00183     if(abs((*iter)->pdg_id())==23){
00184       spinEffectsZ(*iter,weight);
00185     }
00186     if(abs((*iter)->pdg_id())==15){
00187       if(isLastTauinChain(*iter)){
00188         nTaus->Fill(0.5,weight);
00189         int mother  = tauMother(*iter,weight);
00190         int decaychannel = tauDecayChannel(*iter,weight);
00191         tauProngs(*iter, weight);
00192         if(mother>-1){ // exclude B, D and other non-signal decay modes
00193           nPrimeTaus->Fill(0.5,weight);
00194           TauPt->Fill((*iter)->momentum().perp(),weight);
00195           TauEta->Fill((*iter)->momentum().eta(),weight);
00196           TauPhi->Fill((*iter)->momentum().phi(),weight);
00197           rtau(*iter,mother,decaychannel,weight);
00198           photons(*iter,weight);
00199         }
00201         //Adding JAKID and Mass information
00202         //
00203         TauDecay_CMSSW TD;
00204         unsigned int jak_id, TauBitMask;
00205         if(TD.AnalyzeTau((*iter),jak_id,TauBitMask,false,false)){
00206           JAKID->Fill(jak_id,weight);
00207           if(jak_id<=NJAKID){
00208             int tcharge=(*iter)->pdg_id()/abs((*iter)->pdg_id());
00209             std::vector<HepMC::GenParticle*> part=TD.Get_TauDecayProducts();
00210             spinEffects(*iter,mother,jak_id,part,weight);
00211             TLorentzVector LVQ(0,0,0,0);
00212             TLorentzVector LVS12(0,0,0,0);
00213             TLorentzVector LVS13(0,0,0,0);
00214             TLorentzVector LVS23(0,0,0,0);
00215             bool haspart1=false;
00216             for(unsigned int i=0;i<part.size();i++){
00217               if(TD.isTauFinalStateParticle(part.at(i)->pdg_id()) &&
00218                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_e &&
00219                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_mu &&
00220                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_tau ){
00221                 TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00222                 LVQ+=LV;
00223                 if(jak_id==TauDecay::JAK_A1_3PI ||
00224                    jak_id==TauDecay::JAK_KPIK ||
00225                    jak_id==TauDecay::JAK_KPIPI
00226                    ){
00227                   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) ){
00228                     LVS13+=LV;
00229                     LVS23+=LV;
00230                   }
00231                   else{
00232                     LVS12+=LV;
00233                     if(!haspart1 && ((jak_id==TauDecay::JAK_A1_3PI)  || (jak_id!=TauDecay::JAK_A1_3PI && abs(part.at(i)->pdg_id())==PdtPdgMini::K_plus) )){
00234                       LVS13+=LV;
00235                       haspart1=true;
00236                     }
00237                     else{
00238                       LVS23+=LV;
00239                     }
00240                   }
00241                 }
00242               }
00243             }
00244             part.clear();
00245             JAKInvMass.at(jak_id).at(0)->Fill(LVQ.M(),weight);
00246             if(jak_id==TauDecay::JAK_A1_3PI ||
00247                jak_id==TauDecay::JAK_KPIK ||
00248                jak_id==TauDecay::JAK_KPIPI
00249                ){
00250               JAKInvMass.at(jak_id).at(1)->Fill(LVS13.M(),weight);
00251               JAKInvMass.at(jak_id).at(2)->Fill(LVS23.M(),weight);
00252               JAKInvMass.at(jak_id).at(3)->Fill(LVS12.M(),weight);
00253             }
00254           }
00255         }
00256         else{
00257           JAKID->Fill(jak_id,weight);  
00258         }
00259       }
00260     }
00261   }
00262   delete myGenEvent;
00263 }//analyze
00264 
00265 const HepMC::GenParticle* TauValidation::GetMother(const HepMC::GenParticle* tau){
00266   if ( tau->production_vertex() ) {
00267     HepMC::GenVertex::particle_iterator mother;
00268     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00269       if((*mother)->pdg_id() == tau->pdg_id()) return GetMother(*mother);
00270       return (*mother);
00271     }
00272   }
00273   return tau;
00274 }
00275 
00276 
00277 const std::vector<HepMC::GenParticle*> TauValidation::GetMothers(const HepMC::GenParticle* boson){
00278   std::vector<HepMC::GenParticle*> mothers;
00279   if ( boson->production_vertex() ) {
00280     HepMC::GenVertex::particle_iterator mother;
00281     for (mother = boson->production_vertex()->particles_begin(HepMC::parents); mother!= boson->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00282       if((*mother)->pdg_id() == boson->pdg_id()) return GetMothers(*mother);
00283       mothers.push_back(*mother);
00284     }
00285   }
00286   return mothers;
00287 }
00288 
00289 
00290 int TauValidation::findMother(const HepMC::GenParticle* tau){
00291   int mother_pid = 0;
00292   if ( tau->production_vertex() ) {
00293     HepMC::GenVertex::particle_iterator mother;
00294     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00295       mother_pid = (*mother)->pdg_id();
00296       if(mother_pid == tau->pdg_id()) return findMother(*mother); //mother_pid = -1; Make recursive to look for last tau in chain
00297     }
00298   }
00299   return mother_pid;
00300 }
00301 
00302 
00303 bool TauValidation::isLastTauinChain(const HepMC::GenParticle* tau){
00304   if ( tau->end_vertex() ) {
00305     HepMC::GenVertex::particle_iterator dau;
00306     for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
00307       int dau_pid = (*dau)->pdg_id();
00308       if(dau_pid == tau->pdg_id()) return false;
00309     }
00310   }
00311   return true;
00312 }
00313 
00314 
00315 void TauValidation::findTauList(const HepMC::GenParticle* tau,std::vector<const HepMC::GenParticle*> &TauList){
00316   TauList.insert(TauList.begin(),tau);
00317   if ( tau->production_vertex() ) {
00318     HepMC::GenVertex::particle_iterator mother;
00319     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents);mother++) {
00320       if((*mother)->pdg_id() == tau->pdg_id()){
00321         findTauList(*mother,TauList);
00322       }
00323     }
00324   }
00325 }
00326 
00327 void TauValidation::findFSRandBrem(const HepMC::GenParticle* p, bool doBrem, std::vector<const HepMC::GenParticle*> &ListofFSR,
00328                                   std::vector<const HepMC::GenParticle*> &ListofBrem){
00329   // 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. 
00330   if(abs(p->pdg_id())==15){
00331     if(isLastTauinChain(p)){ doBrem=true;}
00332     else{ doBrem=false;}
00333   }
00334     int photo_ID=22;
00335   if ( p->end_vertex() ) {
00336     HepMC::GenVertex::particle_iterator dau;
00337     for (dau = p->end_vertex()->particles_begin(HepMC::children); dau!= p->end_vertex()->particles_end(HepMC::children); dau++ ) {
00338       //if(doBrem) std::cout << "true " << (*dau)->pdg_id() << " "  << std::endl;
00339       //else std::cout << "false " << (*dau)->pdg_id() << " " << std::endl;
00340       if(abs((*dau)->pdg_id()) == abs(photo_ID) && !doBrem){ListofFSR.push_back(*dau);}
00341       if(abs((*dau)->pdg_id()) == abs(photo_ID) && doBrem){ListofBrem.push_back(*dau);}
00342       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*/){
00343         findFSRandBrem(*dau,doBrem,ListofFSR,ListofBrem);
00344       }
00345     }
00346   }
00347 }
00348 
00349 
00350 
00351 void TauValidation::FindPhotosFSR(const HepMC::GenParticle* p,std::vector<const HepMC::GenParticle*> &ListofFSR,double &BosonScale){
00352   BosonScale=0.0;
00353   const HepMC::GenParticle* m=GetMother(p);
00354   double mother_pid=m->pdg_id();
00355   if(m->end_vertex() && mother_pid!=p->pdg_id()){
00356     HepMC::GenVertex::particle_iterator dau;
00357     for (dau = m->end_vertex()->particles_begin(HepMC::children); dau!= m->end_vertex()->particles_end(HepMC::children); dau++ ) {
00358       int dau_pid = (*dau)->pdg_id();
00359       if(fabs(dau_pid) == 22) {
00360         ListofFSR.push_back(*dau);
00361       }
00362     }
00363   }
00364   if(abs(mother_pid) == 24) BosonScale=1.0; // W
00365   if(abs(mother_pid) == 23) BosonScale=2.0; // Z;
00366   if(abs(mother_pid) == 22) BosonScale=2.0; // gamma;
00367   if(abs(mother_pid) == 25) BosonScale=2.0; // HSM;
00368   if(abs(mother_pid) == 35) BosonScale=2.0; // H0;
00369   if(abs(mother_pid) == 36) BosonScale=2.0; // A0;
00370   if(abs(mother_pid) == 37) BosonScale=1.0; //Hpm;
00371 }
00372 
00373 
00374 int TauValidation::tauMother(const HepMC::GenParticle* tau, double weight){
00375 
00376   if(abs(tau->pdg_id()) != 15 ) return -3;
00377   
00378   int mother_pid = findMother(tau);
00379   if(mother_pid == -2) return -2;
00380   
00381   int label = other;
00382   if(abs(mother_pid) == 24) label = W;
00383   if(abs(mother_pid) == 23) label = Z;
00384   if(abs(mother_pid) == 22) label = gamma;
00385   if(abs(mother_pid) == 25) label = HSM;
00386   if(abs(mother_pid) == 35) label = H0;
00387   if(abs(mother_pid) == 36) label = A0;
00388   if(abs(mother_pid) == 37) label = Hpm;
00389   
00390   int mother_shortpid=(abs(mother_pid)%10000);
00391   if(mother_shortpid>500 && mother_shortpid<600 )label = B;
00392   if(mother_shortpid>400 && mother_shortpid<500)label = D;
00393   TauMothers->Fill(label,weight);
00394   if(label==B || label == D || label == other) return -1;
00395   
00396   return mother_pid;
00397 }
00398 
00399 int TauValidation::tauProngs(const HepMC::GenParticle* tau, double weight){
00400   int nProngs = 0;
00401   if ( tau->end_vertex() ) {
00402     HepMC::GenVertex::particle_iterator des;
00403     for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00404         des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00405       int pid = (*des)->pdg_id();
00406       if(abs(pid) == 15) return tauProngs(*des, weight);
00407       if((*des)->status() != 1) continue; // dont count unstable particles
00408       
00409       const HepPDT::ParticleData*  pd = fPDGTable->particle((*des)->pdg_id ());
00410       int charge = (int) pd->charge();
00411       if(charge == 0) continue;
00412       nProngs++;
00413     }
00414   }
00415   TauProngs->Fill(nProngs,weight);
00416   return nProngs;
00417 }
00418 
00419 int TauValidation::tauDecayChannel(const HepMC::GenParticle* tau, double weight){
00420 
00421   int channel = undetermined;
00422   if(tau->status() == 1) channel = stable;
00423   int allCount   = 0,
00424     eCount     = 0,
00425     muCount    = 0,
00426     pi0Count   = 0,
00427     piCount    = 0,
00428     rhoCount   = 0,
00429     a1Count    = 0,
00430     KCount     = 0,
00431     KstarCount = 0;
00432   
00433   if ( tau->end_vertex() ) {
00434     HepMC::GenVertex::particle_iterator des;
00435     for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00436         des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00437       int pid = (*des)->pdg_id();
00438       if(abs(pid) == 15) return tauDecayChannel(*des,weight);
00439       
00440       allCount++;
00441       if(abs(pid) == 11)    eCount++;
00442       if(abs(pid) == 13)    muCount++;
00443       if(abs(pid) == 111)   pi0Count++;
00444       if(abs(pid) == 211)   piCount++;
00445       if(abs(pid) == 213)   rhoCount++;
00446       if(abs(pid) == 20213) a1Count++;
00447       if(abs(pid) == 321)   KCount++;
00448       if(abs(pid) == 323)   KstarCount++;
00449       
00450     }
00451   }
00452   // resonances  
00453   if(KCount >= 1)     channel = K;
00454   if(KstarCount >= 1) channel = Kstar;
00455   if(a1Count >= 1)    channel = a1;
00456   if(rhoCount >= 1)   channel = rho;
00457   if(channel!=undetermined && weight!=0.0) TauDecayChannels->Fill(channel,weight);
00458   
00459   // final state products
00460   if(piCount == 1 && pi0Count == 0) channel = pi;
00461   if(piCount == 1 && pi0Count == 1) channel = pi1pi0;
00462   if(piCount == 1 && pi0Count > 1)  channel = pinpi0;
00463   if(piCount == 3 && pi0Count == 0) channel = tripi;
00464   if(piCount == 3 && pi0Count > 0)  channel = tripinpi0;
00465   if(eCount == 1)                   channel = electron;
00466   if(muCount == 1)                  channel = muon;
00467   if(weight!=0.0) TauDecayChannels->Fill(channel,weight);
00468   return channel;
00469 }
00470 
00471 
00472 void TauValidation::rtau(const HepMC::GenParticle* tau,int mother, int decay, double weight){
00473 
00474         if(decay != pi1pi0) return; // polarization only for 1-prong hadronic taus with one neutral pion to make a clean case
00475 
00476         if(tau->momentum().perp() < tauEtCut) return; // rtau visible only for boosted taus
00477         
00478         double rTau = 0;
00479         double ltrack = leadingPionMomentum(tau, weight);
00480         double visibleTauE = visibleTauEnergy(tau);
00481 
00482         if(visibleTauE != 0) rTau = ltrack/visibleTauE;
00483 
00484         if(abs(mother) == 24) TauRtauW->Fill(rTau,weight);
00485         if(abs(mother) == 37) TauRtauHpm->Fill(rTau,weight); 
00486 }
00487 
00488 void TauValidation::spinEffects(const HepMC::GenParticle* tau,int mother, int decay, std::vector<HepMC::GenParticle*> &part,double weight){
00489   if(decay == TauDecay::JAK_PION || decay == TauDecay::JAK_MUON || decay == TauDecay::JAK_ELECTRON){  // polarization only for 1-prong hadronic taus with no neutral pions
00490     TLorentzVector momP4 = motherP4(tau);
00491     TLorentzVector pionP4 = leadingPionP4(tau);
00492     pionP4.Boost(-1*momP4.BoostVector());
00493     double energy = pionP4.E()/(momP4.M()/2);
00494     if(decay == TauDecay::JAK_PION){
00495       if(abs(mother) == 24) TauSpinEffectsW_X->Fill(energy,weight);     
00496       if(abs(mother) == 37) TauSpinEffectsHpm_X->Fill(energy,weight);
00497     }
00498     if(decay == TauDecay::JAK_MUON){
00499       if(abs(mother) == 24) TauSpinEffectsW_muX->Fill(energy,weight);
00500       if(abs(mother) == 37) TauSpinEffectsHpm_muX->Fill(energy,weight);
00501     }
00502     if(decay == TauDecay::JAK_ELECTRON){
00503       if(abs(mother) == 24) TauSpinEffectsW_eX->Fill(energy,weight);
00504       if(abs(mother) == 37) TauSpinEffectsHpm_eX->Fill(energy,weight);
00505     }
00506 
00507   }
00508   else if(decay==TauDecay::JAK_RHO_PIPI0){
00509     TLorentzVector rho(0,0,0,0),pi(0,0,0,0);
00510     for(unsigned int i=0;i<part.size();i++){
00511       TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00512       if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi_plus){pi+=LV; rho+=LV;}
00513       if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi0){rho+=LV;}
00514     }
00515     if(abs(mother) == 24) TauSpinEffectsW_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
00516     if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
00517   }
00518   else if(decay==TauDecay::JAK_A1_3PI){ // only for pi2pi0 for now
00519     TLorentzVector a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0); 
00520     int nplus(0),nminus(0);
00521     for(unsigned int i=0;i<part.size();i++){
00522       TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00523       if(part.at(i)->pdg_id()==PdtPdgMini::pi_plus){ pi_p+=LV; a1+=LV; nplus++;}
00524       if(part.at(i)->pdg_id()==PdtPdgMini::pi_minus){pi_m+=LV; a1+=LV; nminus++;}
00525     }
00526     double gamma=0;
00527     if(nplus+nminus==3 && nplus==1)  gamma=2*pi_p.Pt()/a1.Pt()-1;
00528     if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.Pt()/a1.Pt()-1;
00529     else{
00530       pi_p+=pi_m; gamma=2*pi_p.Pt()/a1.Pt()-1;
00531     }
00532     if(abs(mother) == 24) TauSpinEffectsW_UpsilonA1->Fill(gamma,weight);
00533     if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonA1->Fill(gamma,weight);
00534   }
00535 }
00536 
00537 void TauValidation::spinEffectsZ(const HepMC::GenParticle* boson, double weight){
00538 
00539   TLorentzVector tautau(0,0,0,0);
00540   TLorentzVector pipi(0,0,0,0);
00541   TLorentzVector taum(0,0,0,0);
00542   int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
00543   double x1(0),x2(0); 
00544   TLorentzVector Zboson(boson->momentum().px(),boson->momentum().py(),boson->momentum().pz(),boson->momentum().e());
00545   if ( boson->end_vertex() ) {
00546     HepMC::GenVertex::particle_iterator des;
00547     for(des = boson->end_vertex()->particles_begin(HepMC::children);    des!= boson->end_vertex()->particles_end(HepMC::children);++des ) {
00548       int pid = (*des)->pdg_id();
00549       if(abs(findMother(*des)) != 15 &&  abs(pid) == 15 && (tauDecayChannel(*des) == pi || tauDecayChannel(*des) == muon || tauDecayChannel(*des) == electron )){
00550         if(tauDecayChannel(*des) == pi)nSinglePionDecays++;
00551         if(tauDecayChannel(*des) == muon)nSingleMuonDecays++;
00552         if(tauDecayChannel(*des) == electron)nSingleElectronDecays++;
00553         TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
00554         tautau += LVtau;
00555         TLorentzVector LVpi=leadingPionP4(*des);
00556         pipi+=LVpi;
00557         const HepPDT::ParticleData*  pd = fPDGTable->particle((*des)->pdg_id ());
00558         int charge = (int) pd->charge();
00559         LVtau.Boost(-1*Zboson.BoostVector());
00560         LVpi.Boost(-1*Zboson.BoostVector());
00561         if(tauDecayChannel(*des) == pi){
00562           if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_X->Fill(LVpi.P()/LVtau.E(),weight);
00563           if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_X->Fill(LVpi.P()/LVtau.E(),weight);
00564         }
00565         if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
00566         else{ x2=LVpi.P()/LVtau.E();}
00567      }
00568     }
00569   }
00570   if(nSingleMuonDecays==2){
00571     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_muX->Fill(x1,weight);
00572     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_muX->Fill(x1,weight);
00573   }
00574   if(nSingleElectronDecays==2){
00575     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_eX->Fill(x1,weight);
00576     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_eX->Fill(x1,weight);
00577   }
00578   if(nSinglePionDecays == 2 && tautau.M()!= 0) {
00579     for(int i=0;i<zsbins;i++){
00580       double zslow=((double)i)*(zsmax-zsmin)/((double)zsbins)+zsmin; 
00581       double zsup=((double)i+1)*(zsmax-zsmin)/((double)zsbins)+zsmin;
00582       double aup=Zstoa(zsup), alow=Zstoa(zslow);
00583       if(x2-x1>alow && x2-x1<aup){
00584         double zs=(zsup+zslow)/2;
00585         if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_Zs->Fill(zs,weight);
00586         if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Zs->Fill(zs,weight);
00587         break;
00588       }
00589     }
00590     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_MVis->Fill(pipi.M()/tautau.M(),weight);
00591     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_MVis->Fill(pipi.M()/tautau.M(),weight);
00592  
00593     if(x1!=0){
00594       const std::vector<HepMC::GenParticle*> m=GetMothers(boson);
00595       int q(0),qbar(0);
00596       TLorentzVector Z(0,0,0,0);
00597       for(unsigned int i=0;i<m.size();i++){
00598         if(m.at(i)->pdg_id()==PdtPdgMini::d      || m.at(i)->pdg_id()==PdtPdgMini::u      ){q++;}
00599         if(m.at(i)->pdg_id()==PdtPdgMini::anti_d || m.at(i)->pdg_id()==PdtPdgMini::anti_u ){qbar++;}
00600       }
00601       if(q==1 && qbar==1){// assume q has largest E (valence vs see quarks) 
00602         if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
00603           if(abs(boson->pdg_id())==PdtPdgMini::Z0)      TauSpinEffectsZ_Xf->Fill(x1,weight);
00604           if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xf->Fill(x1,weight);
00605         }
00606         else{
00607           if(abs(boson->pdg_id())==PdtPdgMini::Z0)      TauSpinEffectsZ_Xb->Fill(x1,weight);
00608           if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xb->Fill(x1,weight);
00609         }
00610       }
00611     }
00612   }
00613 }
00614 
00615 double TauValidation::Zstoa(double zs){
00616   double a=1-sqrt(fabs(1.0-2*fabs(zs)));
00617   if(zs<0){
00618     a*=-1.0;
00619   }
00620   return a;
00621 }
00622 
00623 
00624 double TauValidation::leadingPionMomentum(const HepMC::GenParticle* tau, double weight){
00625         return leadingPionP4(tau).P();
00626 }
00627 
00628 TLorentzVector TauValidation::leadingPionP4(const HepMC::GenParticle* tau){
00629 
00630         TLorentzVector p4(0,0,0,0);
00631 
00632         if ( tau->end_vertex() ) {
00633               HepMC::GenVertex::particle_iterator des;
00634               for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00635                   des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00636                         int pid = (*des)->pdg_id();
00637 
00638                         if(abs(pid) == 15) return leadingPionP4(*des);
00639 
00640                         if(!(abs(pid)==211 || abs(pid)==13 || abs(pid)==11)) continue;
00641  
00642                         if((*des)->momentum().rho() > p4.P()) {
00643                                 p4 = TLorentzVector((*des)->momentum().px(),
00644                                                     (*des)->momentum().py(),
00645                                                     (*des)->momentum().pz(),
00646                                                     (*des)->momentum().e());
00647                         } 
00648                 }
00649         }
00650 
00651         return p4;
00652 }
00653 
00654 TLorentzVector TauValidation::motherP4(const HepMC::GenParticle* tau){
00655   const HepMC::GenParticle* m=GetMother(tau);
00656   return TLorentzVector(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
00657 }
00658 
00659 double TauValidation::visibleTauEnergy(const HepMC::GenParticle* tau){
00660         TLorentzVector p4(tau->momentum().px(),
00661                           tau->momentum().py(),
00662                           tau->momentum().pz(),
00663                           tau->momentum().e());
00664 
00665         if ( tau->end_vertex() ) {
00666               HepMC::GenVertex::particle_iterator des;
00667               for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00668                   des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00669                         int pid = (*des)->pdg_id();
00670 
00671                         if(abs(pid) == 15) return visibleTauEnergy(*des);
00672 
00673                         if(abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
00674                                 p4 -= TLorentzVector((*des)->momentum().px(),
00675                                                      (*des)->momentum().py(),
00676                                                      (*des)->momentum().pz(),
00677                                                      (*des)->momentum().e());
00678                         }
00679                 }
00680         }
00681 
00682         return p4.E();
00683 }
00684 
00685 void TauValidation::photons(const HepMC::GenParticle* tau, double weight){
00686   // Find First tau in chain
00687   std::vector<const HepMC::GenParticle*> TauList;
00688   findTauList(tau,TauList);
00689 
00690   // Get List of Gauge Boson to tau(s) FSR and Brem
00691   bool passedW=false;
00692   std::vector<const HepMC::GenParticle*> ListofFSR;  ListofFSR.clear();
00693   std::vector<const HepMC::GenParticle*> ListofBrem; ListofBrem.clear();
00694   std::vector<const HepMC::GenParticle*> FSR_photos; FSR_photos.clear();
00695   double BosonScale(1);
00696   if(TauList.size()>0){
00697     TauValidation::findFSRandBrem(TauList.at(0),passedW,ListofFSR,ListofBrem);
00698     TauValidation::FindPhotosFSR(TauList.at(0),FSR_photos,BosonScale);
00699 
00700     // Add the Tau Brem. information
00701     TauBremPhotonsN->Fill(ListofBrem.size(),weight);
00702     double photonPtSum=0;
00703     for(unsigned int i=0;i<ListofBrem.size();i++){
00704       photonPtSum+=ListofBrem.at(i)->momentum().perp();
00705       TauBremPhotonsPt->Fill(ListofBrem.at(i)->momentum().perp(),weight);
00706     }
00707     TauBremPhotonsPtSum->Fill(photonPtSum,weight);
00708         
00709     // Now add the Gauge Boson FSR information
00710     if(BosonScale!=0){
00711       TauFSRPhotonsN->Fill(ListofFSR.size(),weight);
00712       photonPtSum=0;
00713       for(unsigned int i=0;i<ListofFSR.size();i++){
00714         photonPtSum+=ListofFSR.at(i)->momentum().perp();
00715         TauFSRPhotonsPt->Fill(ListofFSR.at(i)->momentum().perp(),weight);
00716       }
00717       double FSR_photosSum(0);
00718       for(unsigned int i=0;i<FSR_photos.size();i++){
00719         FSR_photosSum+=FSR_photos.at(i)->momentum().perp();
00720         TauFSRPhotonsPt->Fill(FSR_photos.at(i)->momentum().perp()/BosonScale,weight*BosonScale);
00721       }
00722       TauFSRPhotonsPtSum->Fill(photonPtSum+FSR_photosSum/BosonScale,weight);
00723     }
00724   }
00725 }
00726
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/Validation/EventGenerator/plugins/TauValidation.cc
