lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/tautau/anal
diff -N TauDecayAnal.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TauDecayAnal.java 25 Mar 2009 18:37:16 -0000 1.1
@@ -0,0 +1,497 @@
+package org.lcsim.contrib.Cassell.tautau.anal;
+import org.lcsim.contrib.Cassell.tautau.recon.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.recon.cluster.util.QPhotonClusterEnergyCalculator;
+import org.lcsim.recon.cluster.util.QNeutralHadronClusterEnergyCalculator;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.geometry.IDDecoder;
+import hep.physics.vec.*;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.*;
+import java.util.*;
+import java.io.IOException;
+public class TauDecayAnal extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ int ievt;
+ int ds;
+ String[] dsn = {"+80e-","-80e-"};
+ String[] dsst;
+ String[] idrupst;
+ String[] pmst;
+ int nmax = 1000000;
+ int tauppdg = -15;
+ int taumpdg = 15;
+ String[] chs = {"+","-"};
+ String[] dmodes = {"enunu","mununu","pinu","rhonu","a1nu3ch","Wnu","knu","k*nu","a1nu2pi0","a1nuother",
+ "gammaenunu","gammamununu","gammapinu","gammarhonu","gammaa1nu3ch",
+ "gammaWnu","gammaknu","gammak*nu","gammaa1nu2pi0","gammaa1nuother","unknown"};
+ String[] rmodes = {"elleptonic","muleptonic","1pi+-or1k+","---rho----","a1-2pi01pi","--a1-3pi+-","--other---"};
+ int[] dpdg = {11,13,211,213,20213,20213,15};
+ int[] dtor = {0,1,2,3,5,6,2,6,4,6,0,1,6,6,6,6,6,6,6,6,6};
+ double[][] qcut;
+ double[] qll = {1,3,4};
+ MCParticle[] taumc;
+ TauPairTruthMode ttm;
+ public TauDecayAnal()
+ {
+ ievt = 0;
+ ds = 0;
+ dsst = new String[2];
+ dsst[1] = "both";
+ idrupst = new String[2];
+ idrupst[1] = "both";
+ pmst = new String[2];
+ pmst[1] = "both";
+ ttm = new TauPairTruthMode();
+ add(ttm);
+ taumc = new MCParticle[2];
+ qcut = new double[7][3];
+ qcut[2][1] = .6;
+ qcut[2][0] = .9;
+ qcut[3][2] = .75;
+ qcut[3][1] = .85;
+ qcut[3][0] = .95;
+ qcut[4][2] = .6;
+ qcut[4][1] = .7;
+ qcut[4][0] = .8;
+ qcut[5][0] = .9;
+ }
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+ double wt = (double) event.getWeight();
+ int idrup = event.getIntegerParameters().get("idrup")[0];
+ dsst[0] = dsn[ds];
+ idrupst[0] = "id"+idrup;
+
+
+//
+// Get the truth information
+//
+ int[] gcl = new int[2];
+ int[] gm = new int[2];
+ gcl[0] = ttm.getTauPlusDecayClassInt();
+ gcl[1] = ttm.getTauMinusDecayClassInt();
+ gm[0] = ttm.getTauPlusDecayModeInt();
+ gm[1] = ttm.getTauMinusDecayModeInt();
+ taumc[0] = ttm.getTauPlusMCP();
+ taumc[1] = ttm.getTauMinusMCP();
+ List<ReconstructedParticle> tl = event.get(ReconstructedParticle.class,"TauParticles");
+ ReconstructedParticle tp = null;
+ ReconstructedParticle tm = null;
+ for(ReconstructedParticle p:tl)
+ {
+ if(p.getParticleIDUsed().getPDG() == tauppdg)tp = p;
+ else if(p.getParticleIDUsed().getPDG() == taumpdg)tm = p;
+ }
+ double[] qual = new double[2];
+ qual[0] = tp.getGoodnessOfPID();
+ qual[1] = tm.getGoodnessOfPID();
+ int mp = getMode(tp);
+ int mm = getMode(tm);
+ int[] rm = new int[2];
+ rm[0] = mp;
+ rm[1] = mm;
+ int lp = 0;
+ int lm = 0;
+
+ for(int i=0;i<3;i++)
+ {
+ if(qual[0] < qcut[mp][i])lp++;
+ if(qual[1] < qcut[mm][i])lm++;
+ }
+ int[] ql = new int[2];
+ ql[0] = lp;
+ ql[1] = lm;
+ ReconstructedParticle[] tau = new ReconstructedParticle[2];
+ tau[0] = tp;
+ tau[1] = tm;
+ for(int i=0;i<2;i++)
+ {
+ if(rm[i] != 3)continue;
+ ReconstructedParticle rd = tau[i].getParticles().get(0);
+ double rM = rd.getMass();
+ double rE = rd.getEnergy();
+ double rctst = getcthstRecon(tau[i]);
+ MCParticle tmc = taumc[i];
+ MCParticle d = null;
+ for(MCParticle pp:tmc.getDaughters())
+ {
+ if(pp.getCharge() != 0.)d = pp;
+ }
+ double gM = d.getMass();
+ double gE = d.getEnergy();
+ double gctst = getcthstCheatH(tmc);
+ double gY = 0.;
+ double rY = 0.;
+ double gctstst = 0.;
+ double rctstst = 0.;
+ double rpi0E = 0.;
+ double rpi0M = 0.;
+ ReconstructedParticle rpi0 = null;
+ double gpi0E = 0.;
+ double gpi0M = 0.;
+ MCParticle gpi0 = null;
+ double rcpsi = 0.;
+ double romega = 0.;
+ double gcpsi = 0.;
+ double gomega = 0.;
+ if(rm[i] == 3)
+ {
+ double de = 0.;
+ for(ReconstructedParticle dd:rd.getParticles())
+ {
+ if(de == 0.)de = dd.getEnergy();
+ else de -= dd.getEnergy();
+ if(dd.getCharge() == 0.)rpi0 = dd;
+ }
+ rY = Math.abs(de)/250.;
+ rctstst = getcthststRecon(rd);
+ rcpsi = getcpsiRecon(tau[i]);
+ romega = getomegaRecon(tau[i]);
+ rpi0E = rpi0.getEnergy();
+ rpi0M = rpi0.getMass();
+ }
+ if(gcl[i] == 3)
+ {
+ double de = 0.;
+ for(MCParticle dd:d.getDaughters())
+ {
+ if(de == 0.)de = dd.getEnergy();
+ else de -= dd.getEnergy();
+ if(dd.getCharge() == 0.)gpi0 = dd;
+ }
+ gY = Math.abs(de)/250.;
+ gctstst = getcthststTrue(d);
+ gcpsi = getcpsiCheatH(tmc);
+ gomega = getomegaCheatH(tmc);
+ gpi0E = gpi0.getEnergy();
+ gpi0M = gpi0.getMass();
+// if(gpi0M < .1)throw new RuntimeException("oops");
+ }
+ pmst[0] = "tau"+chs[i];
+ String prem = rmodes[rm[i]];
+ String[] gst = new String[2];
+ gst[1] = "all";
+ if(rm[i] == gcl[i])gst[0] = "corrID";
+ else gst[0] = "wrngID";
+ for(int s1=0;s1<1;s1++)
+ {
+ for(int s2=1;s2<2;s2++)
+ {
+ String pre0 = dsst[s1]+"/"+idrupst[s2]+"/"+prem+"/";
+// String pre0 = idrupst[s2]+"/"+prem+"/";
+ for(int s3=0;s3<2;s3++)
+ {
+ if(ql[i] <= qll[s3])
+ {
+ String pre1 = "qlev "+qll[s3]+"/";
+ double wt2 = 1.;
+ if(s3 == 2)wt2 = qual[i];
+ for(int s4=0;s4<2;s4++)
+ {
+ String pre2 = gst[s4]+"/";
+ String pre = pre0+pre1+pre2;
+ aida.cloud1D(pre+"gen hadron mass",nmax).fill(gM,wt*wt2);
+ aida.cloud1D(pre+"gen hadron energy",nmax).fill(gE,wt*wt2);
+ aida.cloud1D(pre+"gen hadron cos theta*",nmax).fill(gctst,wt*wt2);
+ aida.cloud1D(pre+"rec hadron mass",nmax).fill(rM,wt*wt2);
+ aida.cloud1D(pre+"rec hadron energy",nmax).fill(rE,wt*wt2);
+ aida.cloud1D(pre+"rec hadron cos theta*",nmax).fill(rctst,wt*wt2);
+ aida.cloud1D(pre+"mass residuals",nmax).fill(rM-gM,wt*wt2);
+ aida.cloud1D(pre+"energy residuals",nmax).fill(rE-gE,wt*wt2);
+ aida.cloud1D(pre+"cos theta* residuals",nmax).fill(rctst-gctst,wt*wt2);
+ if(rm[i] == 3)
+ {
+ aida.cloud1D(pre+"gen rho y",nmax).fill(gY,wt*wt2);
+ aida.cloud1D(pre+"rec rho y",nmax).fill(rY,wt*wt2);
+ aida.cloud1D(pre+"rho y residuals",nmax).fill(rY-gY,wt*wt2);
+ aida.cloud1D(pre+"gen rho cos beta",nmax).fill(gctstst,wt*wt2);
+ aida.cloud1D(pre+"rec rho cos beta",nmax).fill(rctstst,wt*wt2);
+ aida.cloud1D(pre+"rho cos beta residuals",nmax).fill(rctstst-gctstst,wt*wt2);
+ aida.cloud1D(pre+"gen rho cos psi",nmax).fill(gcpsi,wt*wt2);
+ aida.cloud1D(pre+"rec rho cos psi",nmax).fill(rcpsi,wt*wt2);
+ aida.cloud1D(pre+"rho cos psi residuals",nmax).fill(rcpsi-gcpsi,wt*wt2);
+ aida.cloud1D(pre+"gen omega",nmax).fill(gomega,wt*wt2);
+ aida.cloud1D(pre+"rec omega",nmax).fill(romega,wt*wt2);
+ aida.cloud1D(pre+"omega residuals",nmax).fill(romega-gomega,wt*wt2);
+ aida.cloud1D(pre+"gen pi0 energy",nmax).fill(gpi0E,wt*wt2);
+ aida.cloud1D(pre+"rec pi0 energy",nmax).fill(rpi0E,wt*wt2);
+ aida.cloud1D(pre+"gen pi0 mass",nmax).fill(gpi0M,wt*wt2);
+ aida.cloud1D(pre+"rec pi0 mass",nmax).fill(rpi0M,wt*wt2);
+ aida.cloud1D(pre+"pi0 mass residuals",nmax).fill(rpi0M-gpi0M,wt*wt2);
+ aida.cloud1D(pre+"pi0 energy residuals",nmax).fill(rpi0E-gpi0E,wt*wt2);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ ievt++;
+ }
+ public int getMode(ReconstructedParticle tau)
+ {
+ ReconstructedParticle d = tau.getParticles().get(0);
+ int pdg = Math.abs(d.getParticleIDUsed().getPDG());
+ int mode = -1;
+ for(int i=0;i<dpdg.length;i++)
+ {
+ if(pdg == dpdg[i])
+ {
+ mode = i;
+ break;
+ }
+ }
+ if(mode == 4)
+ {
+ int nch = 0;
+ for(ReconstructedParticle dd:d.getParticles())
+ {
+ if(dd.getCharge() != 0.)nch++;
+ }
+ if(nch > 1)mode++;
+ }
+ return mode;
+ }
+ protected void suspend()
+ {
+ ds++;
+ }
+ public double getcthstTrue(MCParticle p)
+ {
+ double Et = p.getEnergy();
+ double Mt = p.getMass();
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ MCParticle d = null;
+ for(MCParticle pp:p.getDaughters())
+ {
+ if(pp.getCharge() != 0.)d = pp;
+ }
+ double Eh = d.getEnergy();
+ double Mh = d.getMass();
+ double Ehst = (Mt*Mt + Mh*Mh)/2./Mt;
+ double phst = Math.sqrt(Ehst*Ehst - Mh*Mh);
+ double ctst = (Mt*Eh - Et*Ehst)/(pt*phst);
+ return ctst;
+ }
+ public double getcthstCheatH(MCParticle p)
+ {
+ double Et = 250.;
+ double Mt = 1.777;
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ MCParticle d = null;
+ for(MCParticle pp:p.getDaughters())
+ {
+ if(pp.getCharge() != 0.)d = pp;
+ }
+ double Eh = d.getEnergy();
+ double Mh = d.getMass();
+ double Ehst = (Mt*Mt + Mh*Mh)/2./Mt;
+ double phst = Math.sqrt(Ehst*Ehst - Mh*Mh);
+ double ctst = (Mt*Eh - Et*Ehst)/(pt*phst);
+ return ctst;
+ }
+ public double getcthstRecon(ReconstructedParticle p)
+ {
+ double Et = 250.;
+ double Mt = 1.777;
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ ReconstructedParticle d = null;
+ for(ReconstructedParticle pp:p.getParticles())
+ {
+ if(pp.getCharge() != 0.)d = pp;
+ }
+ double Eh = d.getEnergy();
+ double Mh = d.getMass();
+ double Ehst = (Mt*Mt + Mh*Mh)/2./Mt;
+ double phst = Math.sqrt(Ehst*Ehst - Mh*Mh);
+ double ctst = (Mt*Eh - Et*Ehst)/(pt*phst);
+ return ctst;
+ }
+ public double getcthststTrue(MCParticle p)
+ {
+ double Erho = p.getEnergy();
+ double Mrho = p.getMass();
+ Hep3Vector pr = p.getMomentum();
+ double prho = pr.magnitude();
+ MCParticle pich = null;
+ MCParticle pi0 = null;
+ for(MCParticle pp:p.getDaughters())
+ {
+ if(pp.getCharge() == 0.)pi0 = pp;
+ else pich = pp;
+ }
+ double Epich = pich.getEnergy();
+ double Mpich = pich.getMass();
+ Hep3Vector pc = pich.getMomentum();
+ double ppar = (pr.x()*pc.x()+pr.y()*pc.y()+pr.z()*pc.z())/pr.magnitude();
+ double Mpi0 = pi0.getMass();
+ double gamma = Erho/Mrho;
+ double beta = prho/Erho;
+ double Epichst = (Mrho*Mrho + Mpich*Mpich - Mpi0*Mpi0)/2./Mrho;
+ double ppichst = Math.sqrt(Epichst*Epichst - Mpich*Mpich);
+ double ctst = (gamma*ppar - gamma*beta*Epich)/ppichst;
+ return ctst;
+ }
+ public double getcthststRecon(ReconstructedParticle p)
+ {
+ double Erho = p.getEnergy();
+ double Mrho = p.getMass();
+ Hep3Vector pr = p.getMomentum();
+ double prho = pr.magnitude();
+ ReconstructedParticle pich = null;
+ ReconstructedParticle pi0 = null;
+ for(ReconstructedParticle pp:p.getParticles())
+ {
+ if(pp.getCharge() == 0.)pi0 = pp;
+ else pich = pp;
+ }
+ double Epich = pich.getEnergy();
+ double Mpich = pich.getMass();
+ Hep3Vector pc = pich.getMomentum();
+ double ppar = (pr.x()*pc.x()+pr.y()*pc.y()+pr.z()*pc.z())/pr.magnitude();
+ double Mpi0 = pi0.getMass();
+ double gamma = Erho/Mrho;
+ double beta = prho/Erho;
+ double Epichst = (Mrho*Mrho + Mpich*Mpich - Mpi0*Mpi0)/2./Mrho;
+ double ppichst = Math.sqrt(Epichst*Epichst - Mpich*Mpich);
+ double ctst = (gamma*ppar - gamma*beta*Epich)/ppichst;
+ return ctst;
+ }
+ public double getcpsiTrue(MCParticle tau)
+ {
+ double Et = tau.getEnergy();
+ double Mt = tau.getMass();
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ MCParticle rho = null;
+ for(MCParticle pp:tau.getDaughters())
+ {
+ if(pp.getCharge() != 0.)rho = pp;
+ }
+ double Erho = rho.getEnergy();
+ double Mrho = rho.getMass();
+ double x = Erho/Et;
+ double qsq = Mrho*Mrho;
+ double mtsq = Mt*Mt;
+ double cospsi=(x*(mtsq+qsq)-2*qsq)/((mtsq-qsq)*Math.sqrt(x*x-(qsq/(Et*Et))));
+ return cospsi;
+ }
+ public double getcpsiCheatH(MCParticle tau)
+ {
+ double Et = 250.;
+ double Mt = 1.777;
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ MCParticle rho = null;
+ for(MCParticle pp:tau.getDaughters())
+ {
+ if(pp.getCharge() != 0.)rho = pp;
+ }
+ double Erho = rho.getEnergy();
+ double Mrho = rho.getMass();
+ double x = Erho/Et;
+ double qsq = Mrho*Mrho;
+ double mtsq = Mt*Mt;
+ double cospsi=(x*(mtsq+qsq)-2*qsq)/((mtsq-qsq)*Math.sqrt(x*x-(qsq/(Et*Et))));
+ return cospsi;
+ }
+ public double getcpsiRecon(ReconstructedParticle tau)
+ {
+ double Et = 250.;
+ double Mt = 1.777;
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ ReconstructedParticle rho = null;
+ for(ReconstructedParticle pp:tau.getParticles())
+ {
+ if(pp.getCharge() != 0.)rho = pp;
+ }
+ double Erho = rho.getEnergy();
+ double Mrho = rho.getMass();
+ double x = Erho/Et;
+ double qsq = Mrho*Mrho;
+ double mtsq = Mt*Mt;
+ double cospsi=(x*(mtsq+qsq)-2*qsq)/((mtsq-qsq)*Math.sqrt(x*x-(qsq/(Et*Et))));
+ return cospsi;
+ }
+ public double getomegaTrue(MCParticle tau)
+ {
+ double Et = tau.getEnergy();
+ double Mt = tau.getMass();
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ MCParticle rho = null;
+ for(MCParticle pp:tau.getDaughters())
+ {
+ if(pp.getCharge() != 0.)rho = pp;
+ }
+ double Erho = rho.getEnergy();
+ double Mrho = rho.getMass();
+ double x = Erho/Et;
+ double qsq = Mrho*Mrho;
+ double mtsq = Mt*Mt;
+ double cospsi=(x*(mtsq+qsq)-2*qsq)/((mtsq-qsq)*Math.sqrt(x*x-(qsq/(Et*Et))));
+ double costheta = getcthstTrue(tau);
+ double cosbeta = getcthststTrue(rho);
+ double o1=(-2.+mtsq/qsq+2*(1+mtsq/qsq)*((3*cospsi-1)/2)*((3*cosbeta*cosbeta-1)/2.))*costheta;
+ double o2=3*Math.sqrt(mtsq/qsq)*((3*cosbeta*cosbeta-1)/2)*2*cospsi*Math.sqrt(1-cospsi*cospsi)*
+ Math.sqrt(1-costheta*costheta);
+ double o3=2+mtsq/qsq-2*(1-mtsq/qsq)*((3*cospsi-1)/2)*((3*cosbeta*cosbeta-1)/2.);
+ double omega=(o1+o2)/o3;
+ return omega;
+ }
+ public double getomegaCheatH(MCParticle tau)
+ {
+ double Et = 250.;
+ double Mt = 1.777;
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ MCParticle rho = null;
+ for(MCParticle pp:tau.getDaughters())
+ {
+ if(pp.getCharge() != 0.)rho = pp;
+ }
+ double Erho = rho.getEnergy();
+ double Mrho = rho.getMass();
+ double x = Erho/Et;
+ double qsq = Mrho*Mrho;
+ double mtsq = Mt*Mt;
+ double cospsi=(x*(mtsq+qsq)-2*qsq)/((mtsq-qsq)*Math.sqrt(x*x-(qsq/(Et*Et))));
+ double costheta = getcthstCheatH(tau);
+ double cosbeta = getcthststTrue(rho);
+ double o1=(-2.+mtsq/qsq+2*(1+mtsq/qsq)*((3*cospsi-1)/2)*((3*cosbeta*cosbeta-1)/2.))*costheta;
+ double o2=3*Math.sqrt(mtsq/qsq)*((3*cosbeta*cosbeta-1)/2)*2*cospsi*Math.sqrt(1-cospsi*cospsi)*
+ Math.sqrt(1-costheta*costheta);
+ double o3=2+mtsq/qsq-2*(1-mtsq/qsq)*((3*cospsi-1)/2)*((3*cosbeta*cosbeta-1)/2.);
+ double omega=(o1+o2)/o3;
+ return omega;
+ }
+ public double getomegaRecon(ReconstructedParticle tau)
+ {
+ double Et = 250.;
+ double Mt = 1.777;
+ double pt = Math.sqrt(Et*Et - Mt*Mt);
+ ReconstructedParticle rho = null;
+ for(ReconstructedParticle pp:tau.getParticles())
+ {
+ if(pp.getCharge() != 0.)rho = pp;
+ }
+ double Erho = rho.getEnergy();
+ double Mrho = rho.getMass();
+ double x = Erho/Et;
+ double qsq = Mrho*Mrho;
+ double mtsq = Mt*Mt;
+ double cospsi=(x*(mtsq+qsq)-2*qsq)/((mtsq-qsq)*Math.sqrt(x*x-(qsq/(Et*Et))));
+ double costheta = getcthstRecon(tau);
+ double cosbeta = getcthststRecon(rho);
+ double o1=(-2.+mtsq/qsq+2*(1+mtsq/qsq)*((3*cospsi-1)/2)*((3*cosbeta*cosbeta-1)/2.))*costheta;
+ double o2=3*Math.sqrt(mtsq/qsq)*((3*cosbeta*cosbeta-1)/2)*2*cospsi*Math.sqrt(1-cospsi*cospsi)*
+ Math.sqrt(1-costheta*costheta);
+ double o3=2+mtsq/qsq-2*(1-mtsq/qsq)*((3*cospsi-1)/2)*((3*cosbeta*cosbeta-1)/2.);
+ double omega=(o1+o2)/o3;
+ return omega;
+ }
+}
