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         DecayLength = dbe->book1D("DecayLength","#tau Decay Length", 100 ,0,20);  DecayLength->setAxisTitle("L_{#tau} (mm)");
00077         LifeTime =  dbe->book1D("LifeTime","#tau LifeTime ", 100 ,0,1000E-15);     LifeTime->setAxisTitle("#tau_{#tau}s)");
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     TauSpinEffectsH_pipiAcoplanarity = dbe->book1D("TauSpinEffectsH_pipiAcoplanarity","H Acoplanarity for #pi^{-}#pi^{+}", 50 ,0,2*TMath::Pi()); TauSpinEffectsH_pipiAcoplanarity->setAxisTitle("Acoplanarity"); 
00098 
00099     TauSpinEffectsH_pipiAcollinearity = dbe->book1D("TauSpinEffectsH_pipiAcollinearity","H Acollinearity for #pi^{-}#pi^{+}", 50 ,0,TMath::Pi()); TauSpinEffectsH_pipiAcollinearity->setAxisTitle("Acollinearity");
00100     TauSpinEffectsH_pipiAcollinearityzoom = dbe->book1D("TauSpinEffectsH_pipiAcollinearityzoom","H Acollinearity for #pi^{-}#pi^{+}", 50 ,3,TMath::Pi()); TauSpinEffectsH_pipiAcollinearityzoom->setAxisTitle("Acollinearity");
00101 
00102     TauSpinEffectsZ_MVis   = dbe->book1D("TauSpinEffectsZMVis","Mass of pi+ pi-", 25 ,0,1.1);       TauSpinEffectsZ_MVis->setAxisTitle("M_{#pi^{+}#pi^{-}}");
00103     TauSpinEffectsH_MVis   = dbe->book1D("TauSpinEffectsHMVis","Mass of pi+ pi-", 25 ,0,1.1);       TauSpinEffectsZ_MVis->setAxisTitle("M_{#pi^{+}#pi^{-}}");
00104 
00105     TauSpinEffectsZ_Zs   = dbe->book1D("TauSpinEffectsZZs","Z_{s}", zsbins ,zsmin,zsmax);        TauSpinEffectsZ_Zs->setAxisTitle("Z_{s}");
00106     TauSpinEffectsH_Zs   = dbe->book1D("TauSpinEffectsHZs","Z_{s}", zsbins ,zsmin,zsmax);        TauSpinEffectsZ_Zs->setAxisTitle("Z_{s}");
00107 
00108     TauSpinEffectsZ_X= dbe->book1D("TauSpinEffectsZX","X of #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_X->setAxisTitle("X");
00109     TauSpinEffectsH_X= dbe->book1D("TauSpinEffectsH_X","X of #tau^{-}", 25 ,0,1.0);           TauSpinEffectsH_X->setAxisTitle("X");
00110 
00111     TauSpinEffectsZ_Xf   = dbe->book1D("TauSpinEffectsZXf","X of forward emitted #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_Xf->setAxisTitle("X_{f}");
00112     TauSpinEffectsH_Xf   = dbe->book1D("TauSpinEffectsHXf","X of forward emitted #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_Xf->setAxisTitle("X_{f}");
00113 
00114     TauSpinEffectsZ_Xb   = dbe->book1D("TauSpinEffectsZXb","X of backward emitted #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_Xb->setAxisTitle("X_{b}");
00115     TauSpinEffectsH_Xb   = dbe->book1D("TauSpinEffectsHXb","X of backward emitted #tau^{-}", 25 ,0,1.0);           TauSpinEffectsZ_Xb->setAxisTitle("X_{b}");
00116 
00117     TauSpinEffectsZ_eX   = dbe->book1D("TauSpinEffectsZeX","e energy in Z rest frame", 50 ,0,1);     TauSpinEffectsZ_eX->setAxisTitle("X");
00118     TauSpinEffectsH_eX = dbe->book1D("TauSpinEffectsHeX","e energy in H rest frame", 50 ,0,1); TauSpinEffectsH_eX->setAxisTitle("X");
00119 
00120     TauSpinEffectsZ_muX   = dbe->book1D("TauSpinEffectsZmuX","mu energy in z rest frame", 50 ,0,1);     TauSpinEffectsZ_muX->setAxisTitle("X");
00121     TauSpinEffectsH_muX = dbe->book1D("TauSpinEffectsHmuX","mu energy in H rest frame", 50 ,0,1); TauSpinEffectsH_muX->setAxisTitle("X");
00122 
00123      TauSpinEffectsH_rhorhoAcoplanarityminus = dbe->book1D("TauSpinEffectsH_rhorhoAcoplanarityminus","#phi^{*-} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}<0)", 32 ,0,2*TMath::Pi());     TauSpinEffectsH_rhorhoAcoplanarityminus->setAxisTitle("#phi^{*-} (Acoplanarity)");
00124      TauSpinEffectsH_rhorhoAcoplanarityplus = dbe->book1D("TauSpinEffectsH_rhorhoAcoplanarityplus","#phi^{*+} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}>0)", 32 ,0,2*TMath::Pi());     TauSpinEffectsH_rhorhoAcoplanarityplus->setAxisTitle("#phi^{*+} (Acoplanarity)");
00125 
00126     TauFSRPhotonsN=dbe->book1D("TauFSRPhotonsN","FSR Photons radiating from/with tau (Gauge Boson)", 5 ,-0.5,4.5);
00127     TauFSRPhotonsN->setAxisTitle("N FSR Photons radiating from/with tau");
00128     TauFSRPhotonsPt=dbe->book1D("TauFSRPhotonsPt","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
00129     TauFSRPhotonsPt->setAxisTitle("P_{t} of FSR Photons radiating from/with tau [per tau]");
00130     TauFSRPhotonsPtSum=dbe->book1D("TauFSRPhotonsPtSum","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
00131     TauFSRPhotonsPtSum->setAxisTitle("P_{t} of FSR Photons radiating from/with tau [per tau]");
00132 
00133     TauBremPhotonsN=dbe->book1D("TauBremPhotonsN","Brem. Photons radiating in tau decay", 5 ,-0.5,4.5);
00134     TauBremPhotonsN->setAxisTitle("N FSR Photons radiating from/with tau");
00135     TauBremPhotonsPt=dbe->book1D("TauBremPhotonsPt","Sum Brem Pt ", 100 ,0,100);
00136     TauFSRPhotonsPt->setAxisTitle("P_{t} of Brem. Photons radiating in tau decay");    
00137     TauBremPhotonsPtSum =dbe->book1D("TauBremPhotonsPtSum","Sum of Brem Pt ", 100 ,0,100);
00138     TauFSRPhotonsPtSum->setAxisTitle("Sum P_{t} of Brem. Photons radiating in tau decay");
00139 
00140     JAKID =dbe->book1D("JAKID","JAK ID",NJAKID+1,-0.5,NJAKID+0.5);
00141     for(unsigned int i=0; i<NJAKID+1;i++){
00142       JAKInvMass.push_back(std::vector<MonitorElement *>());
00143       TString tmp="JAKID";
00144       tmp+=i;
00145       JAKInvMass.at(i).push_back(dbe->book1D("M"+tmp,"M_{"+tmp+"} (GeV)", 80 ,0,2.0));
00146       if(i==TauDecay::JAK_A1_3PI ||
00147          i==TauDecay::JAK_KPIK ||
00148          i==TauDecay::JAK_KPIPI ){
00149         JAKInvMass.at(i).push_back(dbe->book1D("M13"+tmp,"M_{13,"+tmp+"} (GeV)", 80 ,0,2.0));
00150         JAKInvMass.at(i).push_back(dbe->book1D("M23"+tmp,"M_{23,"+tmp+"} (GeV)", 80 ,0,2.0));
00151         JAKInvMass.at(i).push_back(dbe->book1D("M12"+tmp,"M_{12,"+tmp+"} (GeV)", 80 ,0,2.0));
00152       }
00153     }
00154   }
00155   
00156   return;
00157 }
00158 
00159 void TauValidation::endJob(){
00160   return;
00161 }
00162 
00163 void TauValidation::beginRun(const edm::Run& iRun,const edm::EventSetup& iSetup)
00164 {
00166   iSetup.getData( fPDGTable );
00167   return;
00168 }
00169 void TauValidation::endRun(const edm::Run& iRun,const edm::EventSetup& iSetup){return;}
00170 void TauValidation::analyze(const edm::Event& iEvent,const edm::EventSetup& iSetup)
00171 { 
00173   edm::Handle<HepMCProduct> evt;
00174   iEvent.getByLabel(hepmcCollection_, evt);
00175 
00176   //Get EVENT
00177   HepMC::GenEvent *myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
00178 
00179   double weight = _wmanager.weight(iEvent);
00180 
00182   /*
00183   edm::Handle<double> WT;
00184   iEvent.getByLabel(edm::InputTag("TauSpinnerGen","TauSpinnerWT"),WT);
00185   weight = 1.0;
00186   if(*(WT.product())>1e-3 && *(WT.product())<=10.0) weight=(*(WT.product()));
00187   else {weight=1.0;}
00188   */
00190 
00191   // find taus
00192   for(HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); iter++) {
00193      if(abs((*iter)->pdg_id())==PdtPdgMini::Z0 || abs((*iter)->pdg_id())==PdtPdgMini::Higgs0){
00194        spinEffectsZH(*iter,weight);
00195     }
00196     if(abs((*iter)->pdg_id())==15){
00197       if(isLastTauinChain(*iter)){
00198         nTaus->Fill(0.5,weight);
00199         int mother  = tauMother(*iter,weight);
00200         int decaychannel = tauDecayChannel(*iter,weight);
00201         tauProngs(*iter, weight);
00202         if(mother>-1){ // exclude B, D and other non-signal decay modes
00203           nPrimeTaus->Fill(0.5,weight);
00204           TauPt->Fill((*iter)->momentum().perp(),weight);
00205           TauEta->Fill((*iter)->momentum().eta(),weight);
00206           TauPhi->Fill((*iter)->momentum().phi(),weight);
00207           rtau(*iter,mother,decaychannel,weight);
00208           photons(*iter,weight);
00209         }
00211         //Adding JAKID and Mass information
00212         //
00213         TauDecay_CMSSW TD;
00214         unsigned int jak_id, TauBitMask;
00215         if(TD.AnalyzeTau((*iter),jak_id,TauBitMask,false,false)){
00216           JAKID->Fill(jak_id,weight);
00217           if(jak_id<=NJAKID){
00218             int tcharge=(*iter)->pdg_id()/abs((*iter)->pdg_id());
00219             std::vector<HepMC::GenParticle*> part=TD.Get_TauDecayProducts();
00220              spinEffectsWHpm(*iter,mother,jak_id,part,weight);
00221             TLorentzVector LVQ(0,0,0,0);
00222             TLorentzVector LVS12(0,0,0,0);
00223             TLorentzVector LVS13(0,0,0,0);
00224             TLorentzVector LVS23(0,0,0,0);
00225             bool haspart1=false;
00226              TVector3 PV,SV;
00227              bool hasDL(false);
00228             for(unsigned int i=0;i<part.size();i++){
00229                if(abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_tau && TD.isTauFinalStateParticle(part.at(i)->pdg_id()) && !hasDL){
00230                  PV=TVector3((*iter)->production_vertex()->point3d().x(),(*iter)->production_vertex()->point3d().y(),(*iter)->production_vertex()->point3d().z());
00231                  SV=TVector3(part.at(i)->production_vertex()->point3d().x(),part.at(i)->production_vertex()->point3d().y(),part.at(i)->production_vertex()->point3d().z());
00232                  TVector3 DL=SV-PV;
00233                  DecayLength->Fill(DL.Mag(),weight);
00234                  double c(2.99792458E8),Ltau(DL.Mag()/1000),beta((*iter)->momentum().rho()/(*iter)->momentum().m());
00235                  LifeTime->Fill( Ltau/(c*beta),weight);
00236                }
00237 
00238               if(TD.isTauFinalStateParticle(part.at(i)->pdg_id()) &&
00239                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_e &&
00240                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_mu &&
00241                  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_tau ){
00242                 TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00243                 LVQ+=LV;
00244                 if(jak_id==TauDecay::JAK_A1_3PI ||
00245                    jak_id==TauDecay::JAK_KPIK ||
00246                    jak_id==TauDecay::JAK_KPIPI
00247                    ){
00248                   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) ){
00249                     LVS13+=LV;
00250                     LVS23+=LV;
00251                   }
00252                   else{
00253                     LVS12+=LV;
00254                     if(!haspart1 && ((jak_id==TauDecay::JAK_A1_3PI)  || (jak_id!=TauDecay::JAK_A1_3PI && abs(part.at(i)->pdg_id())==PdtPdgMini::K_plus) )){
00255                       LVS13+=LV;
00256                       haspart1=true;
00257                     }
00258                     else{
00259                       LVS23+=LV;
00260                     }
00261                   }
00262                 }
00263               }
00264             }
00265             part.clear();
00266             JAKInvMass.at(jak_id).at(0)->Fill(LVQ.M(),weight);
00267             if(jak_id==TauDecay::JAK_A1_3PI ||
00268                jak_id==TauDecay::JAK_KPIK ||
00269                jak_id==TauDecay::JAK_KPIPI
00270                ){
00271               JAKInvMass.at(jak_id).at(1)->Fill(LVS13.M(),weight);
00272               JAKInvMass.at(jak_id).at(2)->Fill(LVS23.M(),weight);
00273               JAKInvMass.at(jak_id).at(3)->Fill(LVS12.M(),weight);
00274             }
00275           }
00276         }
00277         else{
00278           JAKID->Fill(jak_id,weight);  
00279         }
00280       }
00281     }
00282   }
00283   delete myGenEvent;
00284 }//analyze
00285 
00286 const HepMC::GenParticle* TauValidation::GetMother(const HepMC::GenParticle* tau){
00287   if ( tau->production_vertex() ) {
00288     HepMC::GenVertex::particle_iterator mother;
00289     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00290       if((*mother)->pdg_id() == tau->pdg_id()) return GetMother(*mother);
00291       return (*mother);
00292     }
00293   }
00294   return tau;
00295 }
00296 
00297 
00298 const std::vector<HepMC::GenParticle*> TauValidation::GetMothers(const HepMC::GenParticle* boson){
00299   std::vector<HepMC::GenParticle*> mothers;
00300   if ( boson->production_vertex() ) {
00301     HepMC::GenVertex::particle_iterator mother;
00302     for (mother = boson->production_vertex()->particles_begin(HepMC::parents); mother!= boson->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00303       if((*mother)->pdg_id() == boson->pdg_id()) return GetMothers(*mother);
00304       mothers.push_back(*mother);
00305     }
00306   }
00307   return mothers;
00308 }
00309 
00310 
00311 int TauValidation::findMother(const HepMC::GenParticle* tau){
00312   int mother_pid = 0;
00313   if ( tau->production_vertex() ) {
00314     HepMC::GenVertex::particle_iterator mother;
00315     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
00316       mother_pid = (*mother)->pdg_id();
00317       if(mother_pid == tau->pdg_id()) return findMother(*mother); //mother_pid = -1; Make recursive to look for last tau in chain
00318     }
00319   }
00320   return mother_pid;
00321 }
00322 
00323 
00324 bool TauValidation::isLastTauinChain(const HepMC::GenParticle* tau){
00325   if ( tau->end_vertex() ) {
00326     HepMC::GenVertex::particle_iterator dau;
00327     for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
00328       int dau_pid = (*dau)->pdg_id();
00329       if(dau_pid == tau->pdg_id()) return false;
00330     }
00331   }
00332   return true;
00333 }
00334 
00335 
00336 void TauValidation::findTauList(const HepMC::GenParticle* tau,std::vector<const HepMC::GenParticle*> &TauList){
00337   TauList.insert(TauList.begin(),tau);
00338   if ( tau->production_vertex() ) {
00339     HepMC::GenVertex::particle_iterator mother;
00340     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents);mother++) {
00341       if((*mother)->pdg_id() == tau->pdg_id()){
00342         findTauList(*mother,TauList);
00343       }
00344     }
00345   }
00346 }
00347 
00348 void TauValidation::findFSRandBrem(const HepMC::GenParticle* p, bool doBrem, std::vector<const HepMC::GenParticle*> &ListofFSR,
00349                                   std::vector<const HepMC::GenParticle*> &ListofBrem){
00350   // 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. 
00351   if(abs(p->pdg_id())==15){
00352     if(isLastTauinChain(p)){ doBrem=true;}
00353     else{ doBrem=false;}
00354   }
00355     int photo_ID=22;
00356   if ( p->end_vertex() ) {
00357     HepMC::GenVertex::particle_iterator dau;
00358     for (dau = p->end_vertex()->particles_begin(HepMC::children); dau!= p->end_vertex()->particles_end(HepMC::children); dau++ ) {
00359       //if(doBrem) std::cout << "true " << (*dau)->pdg_id() << " "  << std::endl;
00360       //else std::cout << "false " << (*dau)->pdg_id() << " " << std::endl;
00361       if(abs((*dau)->pdg_id()) == abs(photo_ID) && !doBrem){ListofFSR.push_back(*dau);}
00362       if(abs((*dau)->pdg_id()) == abs(photo_ID) && doBrem){ListofBrem.push_back(*dau);}
00363       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*/){
00364         findFSRandBrem(*dau,doBrem,ListofFSR,ListofBrem);
00365       }
00366     }
00367   }
00368 }
00369 
00370 
00371 
00372 void TauValidation::FindPhotosFSR(const HepMC::GenParticle* p,std::vector<const HepMC::GenParticle*> &ListofFSR,double &BosonScale){
00373   BosonScale=0.0;
00374   const HepMC::GenParticle* m=GetMother(p);
00375   double mother_pid=m->pdg_id();
00376   if(m->end_vertex() && mother_pid!=p->pdg_id()){
00377     HepMC::GenVertex::particle_iterator dau;
00378     for (dau = m->end_vertex()->particles_begin(HepMC::children); dau!= m->end_vertex()->particles_end(HepMC::children); dau++ ) {
00379       int dau_pid = (*dau)->pdg_id();
00380       if(fabs(dau_pid) == 22) {
00381         ListofFSR.push_back(*dau);
00382       }
00383     }
00384   }
00385   if(abs(mother_pid) == 24) BosonScale=1.0; // W
00386   if(abs(mother_pid) == 23) BosonScale=2.0; // Z;
00387   if(abs(mother_pid) == 22) BosonScale=2.0; // gamma;
00388   if(abs(mother_pid) == 25) BosonScale=2.0; // HSM;
00389   if(abs(mother_pid) == 35) BosonScale=2.0; // H0;
00390   if(abs(mother_pid) == 36) BosonScale=2.0; // A0;
00391   if(abs(mother_pid) == 37) BosonScale=1.0; //Hpm;
00392 }
00393 
00394 
00395 int TauValidation::tauMother(const HepMC::GenParticle* tau, double weight){
00396 
00397   if(abs(tau->pdg_id()) != 15 ) return -3;
00398   
00399   int mother_pid = findMother(tau);
00400   if(mother_pid == -2) return -2;
00401   
00402   int label = other;
00403   if(abs(mother_pid) == 24) label = W;
00404   if(abs(mother_pid) == 23) label = Z;
00405   if(abs(mother_pid) == 22) label = gamma;
00406   if(abs(mother_pid) == 25) label = HSM;
00407   if(abs(mother_pid) == 35) label = H0;
00408   if(abs(mother_pid) == 36) label = A0;
00409   if(abs(mother_pid) == 37) label = Hpm;
00410   
00411   int mother_shortpid=(abs(mother_pid)%10000);
00412   if(mother_shortpid>500 && mother_shortpid<600 )label = B;
00413   if(mother_shortpid>400 && mother_shortpid<500)label = D;
00414   TauMothers->Fill(label,weight);
00415   if(label==B || label == D || label == other) return -1;
00416   
00417   return mother_pid;
00418 }
00419 
00420 int TauValidation::tauProngs(const HepMC::GenParticle* tau, double weight){
00421   int nProngs = 0;
00422   if ( tau->end_vertex() ) {
00423     HepMC::GenVertex::particle_iterator des;
00424     for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00425         des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00426       int pid = (*des)->pdg_id();
00427       if(abs(pid) == 15) return tauProngs(*des, weight);
00428       if((*des)->status() != 1) continue; // dont count unstable particles
00429       
00430       const HepPDT::ParticleData*  pd = fPDGTable->particle((*des)->pdg_id ());
00431       int charge = (int) pd->charge();
00432       if(charge == 0) continue;
00433       nProngs++;
00434     }
00435   }
00436   TauProngs->Fill(nProngs,weight);
00437   return nProngs;
00438 }
00439 
00440 int TauValidation::tauDecayChannel(const HepMC::GenParticle* tau, double weight){
00441 
00442   int channel = undetermined;
00443   if(tau->status() == 1) channel = stable;
00444   int allCount   = 0,
00445     eCount     = 0,
00446     muCount    = 0,
00447     pi0Count   = 0,
00448     piCount    = 0,
00449     rhoCount   = 0,
00450     a1Count    = 0,
00451     KCount     = 0,
00452     KstarCount = 0;
00453   
00454   if ( tau->end_vertex() ) {
00455     HepMC::GenVertex::particle_iterator des;
00456     for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00457         des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00458       int pid = (*des)->pdg_id();
00459       if(abs(pid) == 15) return tauDecayChannel(*des,weight);
00460       
00461       allCount++;
00462       if(abs(pid) == 11)    eCount++;
00463       if(abs(pid) == 13)    muCount++;
00464       if(abs(pid) == 111)   pi0Count++;
00465       if(abs(pid) == 211)   piCount++;
00466       if(abs(pid) == 213)   rhoCount++;
00467       if(abs(pid) == 20213) a1Count++;
00468       if(abs(pid) == 321)   KCount++;
00469       if(abs(pid) == 323)   KstarCount++;
00470       
00471     }
00472   }
00473   // resonances  
00474   if(KCount >= 1)     channel = K;
00475   if(KstarCount >= 1) channel = Kstar;
00476   if(a1Count >= 1)    channel = a1;
00477   if(rhoCount >= 1)   channel = rho;
00478   if(channel!=undetermined && weight!=0.0) TauDecayChannels->Fill(channel,weight);
00479   
00480   // final state products
00481   if(piCount == 1 && pi0Count == 0) channel = pi;
00482   if(piCount == 1 && pi0Count == 1) channel = pi1pi0;
00483   if(piCount == 1 && pi0Count > 1)  channel = pinpi0;
00484   if(piCount == 3 && pi0Count == 0) channel = tripi;
00485   if(piCount == 3 && pi0Count > 0)  channel = tripinpi0;
00486   if(eCount == 1)                   channel = electron;
00487   if(muCount == 1)                  channel = muon;
00488   if(weight!=0.0) TauDecayChannels->Fill(channel,weight);
00489   return channel;
00490 }
00491 
00492 
00493 void TauValidation::rtau(const HepMC::GenParticle* tau,int mother, int decay, double weight){
00494 
00495         if(decay != pi1pi0) return; // polarization only for 1-prong hadronic taus with one neutral pion to make a clean case
00496 
00497         if(tau->momentum().perp() < tauEtCut) return; // rtau visible only for boosted taus
00498         
00499         double rTau = 0;
00500         double ltrack = leadingPionMomentum(tau, weight);
00501         double visibleTauE = visibleTauEnergy(tau);
00502 
00503         if(visibleTauE != 0) rTau = ltrack/visibleTauE;
00504 
00505         if(abs(mother) == 24) TauRtauW->Fill(rTau,weight);
00506         if(abs(mother) == 37) TauRtauHpm->Fill(rTau,weight); 
00507 }
00508 
00509  void TauValidation::spinEffectsWHpm(const HepMC::GenParticle* tau,int mother, int decay, std::vector<HepMC::GenParticle*> &part,double weight){
00510   if(decay == TauDecay::JAK_PION || decay == TauDecay::JAK_MUON || decay == TauDecay::JAK_ELECTRON){  // polarization only for 1-prong hadronic taus with no neutral pions
00511     TLorentzVector momP4 = motherP4(tau);
00512     TLorentzVector pionP4 = leadingPionP4(tau);
00513     pionP4.Boost(-1*momP4.BoostVector());
00514     double energy = pionP4.E()/(momP4.M()/2);
00515     if(decay == TauDecay::JAK_PION){
00516       if(abs(mother) == 24) TauSpinEffectsW_X->Fill(energy,weight);     
00517       if(abs(mother) == 37) TauSpinEffectsHpm_X->Fill(energy,weight);
00518     }
00519     if(decay == TauDecay::JAK_MUON){
00520       if(abs(mother) == 24) TauSpinEffectsW_muX->Fill(energy,weight);
00521       if(abs(mother) == 37) TauSpinEffectsHpm_muX->Fill(energy,weight);
00522     }
00523     if(decay == TauDecay::JAK_ELECTRON){
00524       if(abs(mother) == 24) TauSpinEffectsW_eX->Fill(energy,weight);
00525       if(abs(mother) == 37) TauSpinEffectsHpm_eX->Fill(energy,weight);
00526     }
00527 
00528   }
00529   else if(decay==TauDecay::JAK_RHO_PIPI0){
00530     TLorentzVector rho(0,0,0,0),pi(0,0,0,0);
00531     for(unsigned int i=0;i<part.size();i++){
00532       TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00533       if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi_plus){pi+=LV; rho+=LV;}
00534       if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi0){rho+=LV;}
00535     }
00536     if(abs(mother) == 24) TauSpinEffectsW_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
00537     if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
00538   }
00539   else if(decay==TauDecay::JAK_A1_3PI){ // only for pi2pi0 for now
00540     TLorentzVector a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0); 
00541     int nplus(0),nminus(0);
00542     for(unsigned int i=0;i<part.size();i++){
00543       TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
00544       if(part.at(i)->pdg_id()==PdtPdgMini::pi_plus){ pi_p+=LV; a1+=LV; nplus++;}
00545       if(part.at(i)->pdg_id()==PdtPdgMini::pi_minus){pi_m+=LV; a1+=LV; nminus++;}
00546     }
00547     double gamma=0;
00548     if(nplus+nminus==3 && nplus==1)  gamma=2*pi_p.Pt()/a1.Pt()-1;
00549     if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.Pt()/a1.Pt()-1;
00550     else{
00551       pi_p+=pi_m; gamma=2*pi_p.Pt()/a1.Pt()-1;
00552     }
00553     if(abs(mother) == 24) TauSpinEffectsW_UpsilonA1->Fill(gamma,weight);
00554     if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonA1->Fill(gamma,weight);
00555   }
00556 }
00557 
00558  void TauValidation::spinEffectsZH(const HepMC::GenParticle* boson, double weight){
00559   TLorentzVector tautau(0,0,0,0);
00560   TLorentzVector pipi(0,0,0,0);
00561   TLorentzVector taum(0,0,0,0);
00562    TLorentzVector taup(0,0,0,0);
00563    TLorentzVector rho_plus,rho_minus,pi_rhominus,pi0_rhominus,pi_rhoplus,pi0_rhoplus,pi_plus,pi_minus;
00564    bool hasrho_minus(false),hasrho_plus(false),haspi_minus(false),haspi_plus(false);
00565   int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
00566   double x1(0),x2(0); 
00567   TLorentzVector Zboson(boson->momentum().px(),boson->momentum().py(),boson->momentum().pz(),boson->momentum().e());
00568   if ( boson->end_vertex() ) {
00569     HepMC::GenVertex::particle_iterator des;
00570     for(des = boson->end_vertex()->particles_begin(HepMC::children);    des!= boson->end_vertex()->particles_end(HepMC::children);++des ) {
00571       int pid = (*des)->pdg_id();
00572       if(abs(findMother(*des)) != 15 &&  abs(pid) == 15 && (tauDecayChannel(*des) == pi || tauDecayChannel(*des) == muon || tauDecayChannel(*des) == electron )){
00573         if(tauDecayChannel(*des) == pi)nSinglePionDecays++;
00574         if(tauDecayChannel(*des) == muon)nSingleMuonDecays++;
00575         if(tauDecayChannel(*des) == electron)nSingleElectronDecays++;
00576         TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
00577         tautau += LVtau;
00578         TLorentzVector LVpi=leadingPionP4(*des);
00579         pipi+=LVpi;
00580         const HepPDT::ParticleData*  pd = fPDGTable->particle((*des)->pdg_id ());
00581         int charge = (int) pd->charge();
00582         LVtau.Boost(-1*Zboson.BoostVector());
00583         LVpi.Boost(-1*Zboson.BoostVector());
00584         if(tauDecayChannel(*des) == pi){
00585           if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_X->Fill(LVpi.P()/LVtau.E(),weight);
00586           if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_X->Fill(LVpi.P()/LVtau.E(),weight);
00587         }
00588         if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
00589         else{ x2=LVpi.P()/LVtau.E();}
00590      }
00591        if(abs(findMother(*des)) != 15 &&  abs(pid) == 15 && (tauDecayChannel(*des) == rho || tauDecayChannel(*des) == pi1pi0) ){
00592          if ( (*des)->end_vertex() ) {
00593            HepMC::GenVertex::particle_iterator tauprod;
00594            TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
00595            if(pid == 15)taum=LVtau; 
00596            if(pid ==-15)taup=LVtau;
00597            for(tauprod = (*des)->end_vertex()->particles_begin(HepMC::descendants); tauprod!= (*des)->end_vertex()->particles_end(HepMC::descendants);++tauprod ) {
00598              int pid_d = (*tauprod)->pdg_id();
00599              if(abs(pid_d)==211 || abs(pid_d)==111){
00600                TLorentzVector LV((*tauprod)->momentum().px(),(*tauprod)->momentum().py(),(*tauprod)->momentum().pz(),(*tauprod)->momentum().e());
00601                if(pid==15){
00602                  hasrho_minus=true;
00603                  if(pid_d==-211 ){ pi_rhominus=LV;}
00604                  if(abs(pid_d)==111 ){ pi0_rhominus=LV;}
00605                }
00606                if(pid==-15){
00607                  hasrho_plus=true;
00608                  if(pid_d==211 ){pi_rhoplus=LV;}
00609                  if(abs(pid_d)==111 ){pi0_rhoplus=LV;} 
00610                }
00611              }
00612            }
00613          }
00614        }
00615        if(abs(findMother(*des)) != 15 &&  abs(pid) == 15 && (tauDecayChannel(*des) == pi) ){
00616          if ( (*des)->end_vertex() ) {
00617            HepMC::GenVertex::particle_iterator tauprod;
00618            TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
00619            if(pid == 15)taum=LVtau;
00620            if(pid ==-15)taup=LVtau;
00621            for(tauprod = (*des)->end_vertex()->particles_begin(HepMC::descendants); tauprod!= (*des)->end_vertex()->particles_end(HepMC::descendants);++tauprod ) {
00622              int pid_d = (*tauprod)->pdg_id();
00623              if(abs(pid_d)==211 || abs(pid_d)==111){
00624                TLorentzVector LV((*tauprod)->momentum().px(),(*tauprod)->momentum().py(),(*tauprod)->momentum().pz(),(*tauprod)->momentum().e());
00625                if(pid==15){
00626                  haspi_minus=true;
00627                  if(pid_d==-211 ){ pi_minus=LV;}
00628                }
00629                if(pid==-15){
00630                  haspi_plus=true;
00631                  if(pid_d==211 ){pi_plus=LV;}
00632                }
00633              }
00634            }
00635          }
00636     }
00637   }
00638    }
00639    if(hasrho_minus && hasrho_plus){
00640      //compute rhorho
00641      rho_minus=pi_rhominus;
00642      rho_minus+=pi0_rhominus;
00643      rho_plus=pi_rhoplus;
00644      rho_plus+=pi0_rhoplus;
00645      TLorentzVector rhorho=rho_minus;rhorho+=rho_plus;
00646 
00647      // boost to rhorho cm
00648      TLorentzVector pi_rhoplusb=pi_rhoplus;     pi_rhoplusb.Boost(-1*rhorho.BoostVector());
00649      TLorentzVector pi0_rhoplusb=pi0_rhoplus;   pi0_rhoplusb.Boost(-1*rhorho.BoostVector());
00650      TLorentzVector pi_rhominusb=pi_rhominus;   pi_rhominusb.Boost(-1*rhorho.BoostVector());
00651      TLorentzVector pi0_rhominusb=pi0_rhominus; pi0_rhominusb.Boost(-1*rhorho.BoostVector());
00652      
00653      // compute n+/-
00654      TVector3 n_plus=pi_rhoplusb.Vect().Cross(pi0_rhoplusb.Vect());
00655      TVector3 n_minus=pi_rhominusb.Vect().Cross(pi0_rhominusb.Vect());
00656 
00657      // compute the acoplanarity
00658      double Acoplanarity=acos(n_plus.Dot(n_minus)/(n_plus.Mag()*n_minus.Mag()));
00659      if(pi_rhominus.Vect().Dot(n_plus)>0){Acoplanarity*=-1;Acoplanarity+=2*TMath::Pi();}
00660 
00661      // now boost to tau frame
00662      pi_rhoplus.Boost(-1*taup.BoostVector());
00663      pi0_rhoplus.Boost(-1*taup.BoostVector());
00664      pi_rhominus.Boost(-1*taum.BoostVector());
00665      pi0_rhominus.Boost(-1*taum.BoostVector());
00666 
00667      // compute y1 and y2
00668      double y1=(pi_rhoplus.E()-pi0_rhoplus.E())/(pi_rhoplus.E()+pi0_rhoplus.E());
00669      double y2=(pi_rhominus.E()-pi0_rhominus.E())/(pi_rhominus.E()+pi0_rhominus.E());
00670 
00671      // fill histograms
00672      if(abs(boson->pdg_id())==PdtPdgMini::Higgs0 && y1*y2<0) TauSpinEffectsH_rhorhoAcoplanarityminus->Fill(Acoplanarity,weight);
00673      if(abs(boson->pdg_id())==PdtPdgMini::Higgs0 && y1*y2>0) TauSpinEffectsH_rhorhoAcoplanarityplus->Fill(Acoplanarity,weight);
00674    }
00675    if(haspi_minus && haspi_plus){
00676      TLorentzVector tauporig=taup;
00677      TLorentzVector taumorig=taum;
00678 
00679      // now boost to Higgs frame
00680      pi_plus.Boost(-1*Zboson.BoostVector());
00681      pi_minus.Boost(-1*Zboson.BoostVector());
00682      
00683      taup.Boost(-1*Zboson.BoostVector());
00684      taum.Boost(-1*Zboson.BoostVector());
00685 
00686      if(abs(boson->pdg_id())==PdtPdgMini::Higgs0){
00687        TauSpinEffectsH_pipiAcollinearity->Fill(acos(pi_plus.Vect().Dot(pi_minus.Vect())/(pi_plus.P()*pi_minus.P())));
00688        TauSpinEffectsH_pipiAcollinearityzoom->Fill(acos(pi_plus.Vect().Dot(pi_minus.Vect())/(pi_plus.P()*pi_minus.P())));
00689      }
00690 
00691      double proj_m=taum.Vect().Dot(pi_minus.Vect())/(taum.P()*taup.P());
00692      double proj_p=taup.Vect().Dot(pi_plus.Vect())/(taup.P()*taup.P());
00693      TVector3 Tau_m=taum.Vect();
00694      TVector3 Tau_p=taup.Vect();
00695      Tau_m*=proj_m;
00696      Tau_p*=proj_p;
00697      TVector3 Pit_m=pi_minus.Vect()-Tau_m;
00698      TVector3 Pit_p=pi_plus.Vect()-Tau_p;
00699 
00700      double Acoplanarity=acos(Pit_m.Dot(Pit_p)/(Pit_p.Mag()*Pit_m.Mag()));
00701      TVector3 n=Pit_p.Cross(Pit_m);
00702      if(n.Dot(Tau_m)/Tau_m.Mag()>0){Acoplanarity*=-1; Acoplanarity+=2*TMath::Pi();}
00703      // fill histograms
00704      if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_pipiAcoplanarity->Fill(Acoplanarity,weight);
00705      taup=tauporig;
00706      taum=taumorig;
00707 
00708    }
00709 
00710 
00711   if(nSingleMuonDecays==2){
00712     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_muX->Fill(x1,weight);
00713     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_muX->Fill(x1,weight);
00714   }
00715   if(nSingleElectronDecays==2){
00716     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_eX->Fill(x1,weight);
00717     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_eX->Fill(x1,weight);
00718   }
00719   if(nSinglePionDecays == 2 && tautau.M()!= 0) {
00720     for(int i=0;i<zsbins;i++){
00721       double zslow=((double)i)*(zsmax-zsmin)/((double)zsbins)+zsmin; 
00722       double zsup=((double)i+1)*(zsmax-zsmin)/((double)zsbins)+zsmin;
00723       double aup=Zstoa(zsup), alow=Zstoa(zslow);
00724       if(x2-x1>alow && x2-x1<aup){
00725         double zs=(zsup+zslow)/2;
00726         if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_Zs->Fill(zs,weight);
00727         if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Zs->Fill(zs,weight);
00728         break;
00729       }
00730     }
00731     if(abs(boson->pdg_id())==PdtPdgMini::Z0)     TauSpinEffectsZ_MVis->Fill(pipi.M()/tautau.M(),weight);
00732     if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_MVis->Fill(pipi.M()/tautau.M(),weight);
00733  
00734     if(x1!=0){
00735       const std::vector<HepMC::GenParticle*> m=GetMothers(boson);
00736       int q(0),qbar(0);
00737       TLorentzVector Z(0,0,0,0);
00738       for(unsigned int i=0;i<m.size();i++){
00739         if(m.at(i)->pdg_id()==PdtPdgMini::d      || m.at(i)->pdg_id()==PdtPdgMini::u      ){q++;}
00740         if(m.at(i)->pdg_id()==PdtPdgMini::anti_d || m.at(i)->pdg_id()==PdtPdgMini::anti_u ){qbar++;}
00741       }
00742       if(q==1 && qbar==1){// assume q has largest E (valence vs see quarks) 
00743         if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
00744           if(abs(boson->pdg_id())==PdtPdgMini::Z0)      TauSpinEffectsZ_Xf->Fill(x1,weight);
00745           if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xf->Fill(x1,weight);
00746         }
00747         else{
00748           if(abs(boson->pdg_id())==PdtPdgMini::Z0)      TauSpinEffectsZ_Xb->Fill(x1,weight);
00749           if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xb->Fill(x1,weight);
00750         }
00751       }
00752     }
00753   }
00754 }
00755 
00756 double TauValidation::Zstoa(double zs){
00757   double a=1-sqrt(fabs(1.0-2*fabs(zs)));
00758   if(zs<0){
00759     a*=-1.0;
00760   }
00761   return a;
00762 }
00763 
00764 
00765 double TauValidation::leadingPionMomentum(const HepMC::GenParticle* tau, double weight){
00766         return leadingPionP4(tau).P();
00767 }
00768 
00769 TLorentzVector TauValidation::leadingPionP4(const HepMC::GenParticle* tau){
00770 
00771         TLorentzVector p4(0,0,0,0);
00772 
00773         if ( tau->end_vertex() ) {
00774               HepMC::GenVertex::particle_iterator des;
00775               for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00776                   des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00777                         int pid = (*des)->pdg_id();
00778 
00779                         if(abs(pid) == 15) return leadingPionP4(*des);
00780 
00781                         if(!(abs(pid)==211 || abs(pid)==13 || abs(pid)==11)) continue;
00782  
00783                         if((*des)->momentum().rho() > p4.P()) {
00784                                 p4 = TLorentzVector((*des)->momentum().px(),
00785                                                     (*des)->momentum().py(),
00786                                                     (*des)->momentum().pz(),
00787                                                     (*des)->momentum().e());
00788                         } 
00789                 }
00790         }
00791 
00792         return p4;
00793 }
00794 
00795 TLorentzVector TauValidation::motherP4(const HepMC::GenParticle* tau){
00796   const HepMC::GenParticle* m=GetMother(tau);
00797   return TLorentzVector(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
00798 }
00799 
00800 double TauValidation::visibleTauEnergy(const HepMC::GenParticle* tau){
00801         TLorentzVector p4(tau->momentum().px(),
00802                           tau->momentum().py(),
00803                           tau->momentum().pz(),
00804                           tau->momentum().e());
00805 
00806         if ( tau->end_vertex() ) {
00807               HepMC::GenVertex::particle_iterator des;
00808               for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
00809                   des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
00810                         int pid = (*des)->pdg_id();
00811 
00812                         if(abs(pid) == 15) return visibleTauEnergy(*des);
00813 
00814                         if(abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
00815                                 p4 -= TLorentzVector((*des)->momentum().px(),
00816                                                      (*des)->momentum().py(),
00817                                                      (*des)->momentum().pz(),
00818                                                      (*des)->momentum().e());
00819                         }
00820                 }
00821         }
00822 
00823         return p4.E();
00824 }
00825 
00826 void TauValidation::photons(const HepMC::GenParticle* tau, double weight){
00827   // Find First tau in chain
00828   std::vector<const HepMC::GenParticle*> TauList;
00829   findTauList(tau,TauList);
00830 
00831   // Get List of Gauge Boson to tau(s) FSR and Brem
00832   bool passedW=false;
00833   std::vector<const HepMC::GenParticle*> ListofFSR;  ListofFSR.clear();
00834   std::vector<const HepMC::GenParticle*> ListofBrem; ListofBrem.clear();
00835   std::vector<const HepMC::GenParticle*> FSR_photos; FSR_photos.clear();
00836   double BosonScale(1);
00837   if(TauList.size()>0){
00838     TauValidation::findFSRandBrem(TauList.at(0),passedW,ListofFSR,ListofBrem);
00839     TauValidation::FindPhotosFSR(TauList.at(0),FSR_photos,BosonScale);
00840 
00841     // Add the Tau Brem. information
00842     TauBremPhotonsN->Fill(ListofBrem.size(),weight);
00843     double photonPtSum=0;
00844     for(unsigned int i=0;i<ListofBrem.size();i++){
00845       photonPtSum+=ListofBrem.at(i)->momentum().perp();
00846       TauBremPhotonsPt->Fill(ListofBrem.at(i)->momentum().perp(),weight);
00847     }
00848     TauBremPhotonsPtSum->Fill(photonPtSum,weight);
00849         
00850     // Now add the Gauge Boson FSR information
00851     if(BosonScale!=0){
00852       TauFSRPhotonsN->Fill(ListofFSR.size(),weight);
00853       photonPtSum=0;
00854       for(unsigned int i=0;i<ListofFSR.size();i++){
00855         photonPtSum+=ListofFSR.at(i)->momentum().perp();
00856         TauFSRPhotonsPt->Fill(ListofFSR.at(i)->momentum().perp(),weight);
00857       }
00858       double FSR_photosSum(0);
00859       for(unsigned int i=0;i<FSR_photos.size();i++){
00860         FSR_photosSum+=FSR_photos.at(i)->momentum().perp();
00861         TauFSRPhotonsPt->Fill(FSR_photos.at(i)->momentum().perp()/BosonScale,weight*BosonScale);
00862       }
00863       TauFSRPhotonsPtSum->Fill(photonPtSum+FSR_photosSum/BosonScale,weight);
00864     }
00865   }
00866 }
00867
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/Validation/EventGenerator/plugins/TauValidation.cc
