lcsim/sandbox/NickSinev/Examples
diff -N CovMatrixPlot.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CovMatrixPlot.java 8 Jun 2010 18:02:45 -0000 1.1
@@ -0,0 +1,510 @@
+import java.util.*;
+import java.io.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.constants.Constants;
+import hep.aida.ITree;
+import java.util.Random;
+import hep.physics.vec.*;
+import hep.aida.*;
+import java.text.*;
+
+
+/**
+ *
+ * @author Nick Sinev
+ */
+public class CovMatrixPlot
+{
+ static int ctc=7;
+ static int ifc=24;
+ static int wtpl=0;
+ static boolean shbns = false;
+ static String fname = "Bruce_finebns_24thmsk_nbc.txt";
+ static String reffnam="Bruce_finebns_nomskl_nbc.txt";
+ static double[][][] cmdiag = new double[5][41][51];
+ static double[] ctbins = new double[41];
+ static double[] pbins = new double[51];
+ static double costhet = 0.;
+ static String ctstr="_cost89_mstrkb4";
+ static String titstr="cos(theta)=0.89";
+ static int ctbi=19;
+ static double[] refer = new double[51];
+ static String conds = " missing 1st vtx endcap layer";
+ static String msstr=" ";
+/*
+*/
+ public static void main(String[] args) throws IOException
+ {
+ switch(ifc)
+ {
+ case 0:
+ msstr="msvtb1";
+ conds="missing 1st vtx bar layer";
+ fname="Bruce_finebns_0thmsk_nbc.txt";
+ break;
+ case 1:
+ msstr="msvtb2";
+ conds="missing 2nd vtx bar layer";
+ fname="Bruce_finebns_1stmsk_nbc.txt";
+ break;
+ case 2:
+ msstr="msvtb3";
+ conds="missing 3rd vtx bar layer";
+ fname="Bruce_finebns_2ndmsk_nbc.txt";
+ break;
+ case 3:
+ msstr="msvtb4";
+ conds="missing 4th vtx bar layer";
+ fname="Bruce_finebns_3rdmsk_nbc.txt";
+ break;
+ case 4:
+ msstr="msvtb5";
+ conds="missing 5th vtx bar layer";
+ fname="Bruce_finebns_4thmsk_nbc.txt";
+ break;
+ case 5:
+ msstr="msvtec1";
+ conds="missing 1st vtx ec layer";
+ fname="Bruce_finebns_5thmsk_nbc.txt";
+ break;
+ case 6:
+ msstr="msvtec2";
+ conds="missing 2nd vtx ec layer";
+ fname="Bruce_finebns_6thmsk_nbc.txt";
+ break;
+ case 7:
+ msstr="msvtec3";
+ conds="missing 3rd vtx ec layer";
+ fname="Bruce_finebns_7thmsk_nbc.txt";
+ break;
+ case 8:
+ msstr="msvtec4";
+ conds="missing 4th vtx ec layer";
+ fname="Bruce_finebns_8thmsk_nbc.txt";
+ break;
+ case 9:
+ msstr="mstrkb1";
+ conds="missing 1st trk bar layer";
+ fname="Bruce_finebns_9thmsk_nbc.txt";
+ break;
+ case 10:
+ msstr="mstrkb2";
+ conds="missing 2nd trk bar layer";
+ fname="Bruce_finebns_10thmsk_nbc.txt";
+ break;
+ case 11:
+ msstr="mstrkb3";
+ conds="missing 3rd trk bar layer";
+ fname="Bruce_finebns_11thmsk_nbc.txt";
+ break;
+ case 12:
+ msstr="mstrkb4";
+ conds="missing 4th trk bar layer";
+ fname="Bruce_finebns_12thmsk_nbc.txt";
+ break;
+ case 13:
+ msstr="mstrkb5";
+ conds="missing last trk bar layer";
+ fname="Bruce_finebns_13thmsk_nbc.txt";
+ break;
+ case 14:
+ msstr="mstrkec1";
+ conds="missing first trk ec layer";
+ fname="Bruce_finebns_14thmsk_nbc.txt";
+ break;
+ case 16:
+ msstr="mstrkec2";
+ conds="missing 2nd trk ec layer";
+ fname="Bruce_finebns_16thmsk_nbc.txt";
+ break;
+ case 18:
+ msstr="mstrkec3";
+ conds="missing 3rd trk ec layer";
+ fname="Bruce_finebns_18thmsk_nbc.txt";
+ break;
+ case 20:
+ msstr="mstrkec4";
+ conds="missing last trk EC layer";
+ fname="Bruce_finebns_20thmsk_nbc.txt";
+ break;
+ case 22:
+ msstr="mstrkfwd1";
+ conds="missing first fwd dsk";
+ fname="Bruce_finebns_22ndmsk_nbc.txt";
+ break;
+ case 23:
+ msstr="mstrkfwd2";
+ conds="missing second fwd dsk";
+ fname="Bruce_finebns_23rdmsk_nbc.txt";
+ break;
+ case 24:
+ msstr="mstrkfwd3";
+ conds="missing third fwd dsk";
+ fname="Bruce_finebns_24thmsk_nbc.txt";
+ break;
+ default:
+ System.out.println("undefined value for input file choise!");
+ return;
+ }
+ switch(ctc)
+ {
+ case 0:
+ ctstr="_cost00_"+msstr;
+ titstr="cos(theta)=0.0";
+ ctbi=0;
+ costhet=0.;
+ break;
+ case 1:
+ ctstr="_cost50_"+msstr;
+ titstr="cos(theta)=0.5";
+ ctbi=6;
+ costhet=0.5;
+ break;
+ case 2:
+ ctstr="_cost75_"+msstr;
+ titstr="cos(theta)=0.75";
+ ctbi=12;
+ costhet=0.749;
+ break;
+ case 3:
+ ctstr="_cost80_"+msstr;
+ titstr="cos(theta)=0.8";
+ ctbi=14;
+ costhet=0.8;
+ break;
+ case 4:
+ ctstr="_cost89_"+msstr;
+ titstr="cos(theta)=0.89";
+ ctbi=19;
+ costhet=0.89;
+ break;
+ case 5:
+ ctstr="_cost95_"+msstr;
+ titstr="cos(theta)=0.95";
+ ctbi=26;
+ costhet=0.95;
+ break;
+ case 6:
+ ctstr="_cost97_"+msstr;
+ titstr="cos(theta)=0.972";
+ ctbi=31;
+ costhet=0.972;
+ break;
+ case 7:
+ ctstr="_cost99_"+msstr;
+ titstr="cos(theta)=0.986";
+ ctbi=37;
+ costhet=0.986;
+ break;
+ default:
+ System.out.println("undefined value for cos theta choise!");
+ return;
+ }
+
+ /*
+ * Now we are going to read file, and can continue only if it was successful
+ */
+ double sinth = Math.sqrt(1.-costhet*costhet);
+ for(wtpl=0; wtpl<5; wtpl++)
+ {
+ boolean file_found = readCovMatrixDiag(reffnam);
+ if(!file_found)
+ {
+ System.out.println("Could not find "+reffnam);
+ return;
+ }
+ if((wtpl==0) && (shbns))
+ {
+ for(int i=0; i<41; i++)
+ System.out.println("bin "+i+" cos theta: "+ctbins[i]);
+ }
+ for(int i=0; i<51; i++)
+ {
+ refer[i]=cmdiag[wtpl][ctbi][i];
+ }
+ file_found=readCovMatrixDiag(fname);
+ if(!file_found)
+ {
+ System.out.println("Could not find "+fname);
+ return;
+ }
+ /*
+ * now work with AIDA
+ */
+
+ AIDA aida = AIDA.defaultInstance();
+ ITree tree = aida.tree();
+ IAnalysisFactory af = IAnalysisFactory.create();
+ IHistogramFactory hf = af.createHistogramFactory(tree);
+ IDataPointSetFactory dpsf = af.createDataPointSetFactory(tree);
+
+ /*
+ * The best way to plot trace is to use DataPointSet object.
+ */
+ String reftitp = "d0";
+ if(wtpl==1) reftitp="phi0";
+ if(wtpl==2) reftitp="rel_omega";
+ if(wtpl==3) reftitp="z0";
+ if(wtpl==4) reftitp="Tan_Lam";
+ String yaxlbl = "dD0 (microns) ";
+ if(wtpl==1) yaxlbl="dphi0 ";
+ if(wtpl==2) yaxlbl="dPt/Pt ";
+ if(wtpl==3) yaxlbl="dz0 (microns) ";
+ if(wtpl==4) yaxlbl="dTl ";
+ String reftitb = " resolution vs full momentum, cos theta = ";
+ String reftit = reftitp+reftitb+costhet;
+ IDataPointSet dpsr = dpsf.create("all layers working",reftit,2);
+ IDataPointSet dpsm = dpsf.create(conds,reftitp+reftitb+conds,2);
+
+ double Pt=10.0;
+ double omega = (5.*Constants.fieldConversion)/Pt;
+// double curad=1./omega;
+// System.out.println("Got curvature radius: "+curad+" for Pt= "+Pt+" omega is "+omega);
+ int ntry = 51;
+ int j=0;
+ double maxv=0.;
+ double minv=100000.;
+ double val=0.;
+ for(int tr=0; tr<ntry; tr++)
+ {
+ double fmom = pbins[j];
+ Pt = fmom*sinth;
+ omega = (5.*Constants.fieldConversion)/Pt;
+ dpsr.addPoint();
+ dpsr.point(j).coordinate(0).setValue(pbins[j]);
+ if(wtpl==0) val=10000.*Math.sqrt(refer[j]);
+ if(wtpl==1) val=Math.sqrt(refer[j]);
+ if(wtpl==2) val=0.1*Math.sqrt(refer[j])/omega;
+ if(wtpl==3) val=10000.*Math.sqrt(refer[j]);
+ if(wtpl==4) val=Math.sqrt(refer[j]);
+ dpsr.point(j).coordinate(1).setValue(val);
+ if(val > maxv) maxv=val;
+ if(val < minv) minv=val;
+ dpsm.addPoint();
+ dpsm.point(j).coordinate(0).setValue( pbins[j]);
+ if(wtpl==0) val=10000.*Math.sqrt(cmdiag[0][ctbi][j]) ; // cm to micron conversion
+ if(wtpl==1) val=Math.sqrt(cmdiag[1][ctbi][j]) ;
+ if(wtpl==2) val= 0.1*Math.sqrt(cmdiag[2][ctbi][j] )/omega;
+ if(wtpl==3) val= 10000.*Math.sqrt(cmdiag[3][ctbi][j] );
+ if(wtpl==4) val= Math.sqrt(cmdiag[4][ctbi][j] );
+ dpsm.point(j).coordinate(1).setValue(val );
+// if(j==30) System.out.println("For fmom= "+fmom+" Pt= "+Pt+" omega= "+omega+" cov.matr. el "+cmdiag[2][ctbi][j]+" plot value "+val);
+ if(val > maxv) maxv=val;
+ if(val < minv) minv=val;
+ j++;
+ }
+ /*
+* Now set up plotting
+*/
+ IPlotterFactory pf = af.createPlotterFactory();
+ IPlotter plotter = pf.create("Plot");
+// plotter.createRegions(1,1);
+ IPlotterRegion reg0 = plotter.region(0);
+ reg0.setXLimits(pbins[0],pbins[50]);
+ reg0.setYLimits(minv,maxv);
+// System.out.println("Setting min y value: "+minv+" and max: "+maxv);
+ IPlotterStyle rstyle = reg0.style();
+ ILegendBoxStyle legst = rstyle.legendBoxStyle();
+ legst.setVisible(true);
+ legst.textStyle().setFontSize(14.);
+ legst.textStyle().setBold(true);
+ rstyle.titleStyle().setVisible(false);
+ rstyle.xAxisStyle().labelStyle().setFontSize(20);
+ rstyle.xAxisStyle().labelStyle().setBold(true);
+ rstyle.xAxisStyle().setLabel("Tracks momentum (GeV/c) ");
+ rstyle.xAxisStyle().tickLabelStyle().setFontSize(16);
+ rstyle.xAxisStyle().tickLabelStyle().setBold(true);
+ rstyle.xAxisStyle().setScaling("log");
+ rstyle.yAxisStyle().labelStyle().setFontSize(20);
+ rstyle.yAxisStyle().labelStyle().setBold(true);
+ rstyle.yAxisStyle().setLabel(yaxlbl);
+ rstyle.yAxisStyle().tickLabelStyle().setFontSize(16);
+ rstyle.yAxisStyle().tickLabelStyle().setBold(true);
+ rstyle.yAxisStyle().setScaling("log");
+ rstyle.statisticsBoxStyle().setVisible(false);
+ rstyle.dataStyle().outlineStyle().setVisible(true);
+ rstyle.dataStyle().outlineStyle().setThickness(4);
+ rstyle.dataStyle().outlineStyle().setColor("Blue");
+ IPlotterStyle style =pf.createPlotterStyle();
+ style.dataStyle().markerStyle().setParameter("size","4");
+ style.dataStyle().markerStyle().setParameter("shape","1");
+ style.dataStyle().markerStyle().setParameter("color","blue");
+ rstyle.titleStyle().setVisible(true);
+ style.titleStyle().textStyle().setFontSize(14);
+ style.titleStyle().textStyle().setBold(true);
+ plotter.region(0).plot(dpsr,style);
+ style.dataStyle().outlineStyle().setColor("red");
+ style.dataStyle().markerStyle().setParameter("color","red");
+ plotter.region(0).plot(dpsm,style);
+ plotter.region(0).setTitle(titstr);
+ // plotter.titleStyle().setVisible(true);
+// plotter.titleStyle().textStyle().setFontSize(14);
+// plotter.titleStyle().textStyle().setBold(true);
+
+ plotter.show();
+// plotter.writeToFile("C:\\analysis\\LOI\\mlreff_"+reftitp+ctstr+".jpg");
+ }
+ }
+
+ private static boolean readCovMatrixDiag(String fname) throws IOException
+ {
+ File cache;
+ cache = new File(fname);
+ if(!cache.exists())
+ {
+ String cacheDir = System.getProperty("user.home");
+ File cachedir = new File(cacheDir);
+ if(cachedir == null)
+ {
+ System.out.println("File "+fname+" does not exist, and user home directory is not found!");
+ return false;
+ }
+ else
+ {
+ File home = new File(cachedir,".cache");
+ cache = new File(home,fname);
+ }
+ }
+ if(cache.exists())
+ {
+ System.out.println("Found file: "+cache.getAbsolutePath());
+ FileInputStream fins = new FileInputStream(cache);
+ BufferedReader r = new BufferedReader(new InputStreamReader(fins));
+ String rdstr;
+ boolean indiagrd=false;
+ int ind1=0;
+ int ind2=0;
+ int lnmbr=0;
+ int fldst=0;
+ int fldend=0;
+ String fld = null;
+ int nctbns=0;
+ int npbns=0;
+ double[] dvals= new double[4];
+ int ndv=0;
+ int ctind=0;
+ int pind=0;
+ double ctval=0.;
+ while((rdstr=r.readLine()) != null)
+ {
+ if(rdstr.indexOf("entry") != -1)
+ {
+ fldst = rdstr.indexOf("entry");
+ fldst+=5;
+ fldend=rdstr.indexOf(",");
+ fld = rdstr.substring(fldst,fldend);
+ ind1 = Integer.parseInt(fld.trim());
+ fldst=fldend;
+ fldst++;
+ fld=rdstr.substring(fldst,rdstr.length());
+ ind2 = Integer.parseInt(fld.trim());
+ if(ind1==ind2) indiagrd=true;
+ lnmbr=0;
+ ctind=0;
+ pind=0;
+ rdstr=r.readLine();
+ if(rdstr==null) break;
+ }
+ if(!indiagrd) continue;
+ if((rdstr.trim()).length()<2) continue;
+ lnmbr++;
+ if(lnmbr==1)
+ {
+ nctbns=Integer.parseInt(rdstr.trim());
+ if(nctbns != 41) System.out.println("Number of cos theta bins "+nctbns+" != 41");
+ }
+ if(lnmbr==2)
+ {
+ npbns=Integer.parseInt(rdstr.trim());
+ if(npbns!=51) System.out.println("Number of momentum bins "+npbns+" != 51");
+ }
+ if(lnmbr==3)
+ {
+ ctval=Double.parseDouble(rdstr.trim());
+ ctbins[ctind]=ctval;
+ pind=0;
+ }
+ if(lnmbr>3)
+ {
+ if(pind < npbns)
+ {
+ ndv=getDoublesFromString(rdstr,dvals);
+ if(ndv != 2) System.out.println("number of double values in string "+ndv+" not equal 2!");
+ if(ndv==2)
+ {
+ pbins[pind]=dvals[0];
+ cmdiag[ind1-1][ctind][pind]=dvals[1];
+ pind++;
+ }
+ }
+ if(pind>=npbns)
+ {
+ ctind++;
+ if(ctind>=nctbns) indiagrd=false;
+ lnmbr=2;
+ continue;
+ }
+ } // end if(lnmbr>3)
+ } // end while
+ return true;
+ } // end if(cashe.exists()
+ return false;
+}
+
+ private static int getDoublesFromString(String str, double[] buf)
+ {
+ int blen = buf.length;
+ int nv = 0;
+ int inofb = str.indexOf("}");
+ String clstr =str;
+ if(inofb !=-1) clstr = str.substring(0,inofb);
+ String[] tokens = clstr.split("\\s+");
+ int ntp = tokens.length;
+ for(int i=0; i<ntp; i++)
+ {
+ try
+ {
+ double dval = Double.parseDouble(tokens[i]);
+ if(nv<blen)
+ buf[nv]=dval;
+ nv++;
+ }
+ catch(NumberFormatException e)
+ {
+ }
+ finally
+ {
+ }
+ }
+ return nv;
+ }
+
+ private static int getIntegersFromString(String str, int[] buf)
+ {
+ int blen = buf.length;
+ int nv = 0;
+ int inofb = str.indexOf("}");
+ String clstr =str;
+ if(inofb !=-1) clstr = str.substring(0,inofb);
+ String[] tokens = clstr.split("\\s+");
+ int ntp = tokens.length;
+ for(int i=0; i<ntp; i++)
+ {
+ try
+ {
+ int ival = Integer.parseInt(tokens[i]);
+ if(nv<blen)
+ buf[nv]=ival;
+ nv++;
+ }
+ catch(NumberFormatException e)
+ {
+ }
+ finally
+ {
+ }
+ }
+ return nv;
+ }
+
+}
\ No newline at end of file
