Skip to content

Commit b088fe6

Browse files
raffaelladevitaraffaelladevita
authored
DC alignment updates (#18)
* added beam offset analysis * used moment at R1 * saving new offset histograms to the histogram file * improving offset fit and reverting to use track vertex variables * added theta threshold function with offset * changing function names to be sector dependent * restore time residual plots, update to skimming script * restore time residual plots, attempt2 * restoring writing of time residuals to file * reading time histograms when reanalyzing or refitting * restoring vertex constraint * set a 10 count threshold on the scattering chamber exit window peak for including the fitted position in the global fit * improved vertex fit for new cryotarget * added time residuals calibration plots, switched to newer coatjava version, added target length constraint, changed global transformation to translation only * modified reading and writing of histogram file to save fitted functions * added fit of integrated vertex distribution, setting of residual stat box accuracy * renamed histograms * disabled jminuit log during histogram fitting * added checks on histograms fit convergence and code cleanups * restored log file * now selecting status=0 hits, added wire plots * r1 translations are now global, changing theta to traj point * bug fix * added option to initialize misalignment fit from CCDB table (-init), previous iteration table option changed to -previous * added command-line options to choose frame and global/loval R1 translation * cleanup * increased R3 residual error by x2 and added alpha plots * changed alpha to reducedalpha * fix to global transformation handling and more time plots * more plots and improved fitting * improved DC skimming scripts * restoring time residual bins for angular bins * adding doca cuts * improving hit selection * added alpha cuts on hits and double gaussian fit of residuals * improving hit selection again * redo failed fits * added command line options to turn hit cuts on/off * fix setting of doca cuts and added more logging * added command line option to set the ecal cut, plotting pulls instead of residuals in before/after, better labels and printouts * added vertex fit for Spring18 * fixing the target length in the RG-A Spring18 fits * changing moller cone distance to FTon * added angular plot of time residuals * do not refit time residual if fit is stored * creating timeValues array even if tres is false * version bump * added time residuals for left right hits * making LR time residual plots layer dependent * saving LR time residuals to the histogram file * saving LR time residuals to the histogram file * added region-by-region before/after plots * file name now shown in the jframe title and code cleanups * more vertex plots in kinematics.groovy * added RG-K vertex fit function, excluding measurements with failed fits, storing failed fit information, analyze now saves the histograms, using inputlist for analyze and fit mode * now really switching -fit parser options to use inputlist * updated .gitignore with dc specific files * plot only non empty graphs when running only with r0 * fixed issue with skim script * version bump
1 parent eb38e3d commit b088fe6

File tree

7 files changed

+263
-50
lines changed

7 files changed

+263
-50
lines changed

‎.gitignore‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,13 @@ target/
33
*.class
44
*.jar
55

6+
# DC-specific files.
7+
dc/java/dc-alignment/*.hipo
8+
dc/java/dc-alignment/dc-alignment.log*
9+
dc/java/dc-alignment/*.sqlite
10+
dc/java/dc-alignment/*.txt
11+
dc/utilities/*.hipo
12+
613
# FMT-specific files.
714
fmt/*.hipo
815
fmt/error_report.txt

‎dc/java/dc-alignment/pom.xml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
<modelVersion>4.0.0</modelVersion>
44
<groupId>org.clas.detector</groupId>
55
<artifactId>dc-alignment</artifactId>
6-
<version>2.1</version>
6+
<version>2.2</version>
77
<packaging>jar</packaging>
88

99
<properties>

‎dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Alignment.java‎

Lines changed: 98 additions & 32 deletions
Large diffs are not rendered by default.

‎dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Constants.java‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,7 @@ public class Constants {
4242
public static int VDFBINS = 200;
4343
public static double VDFMIN = -5.0;
4444
public static double VDFMAX = 5.0;
45+
public static double CHI2MAX = 2.5;
4546

4647

4748
// global fit
@@ -74,6 +75,7 @@ public class Constants {
7475

7576

7677
// measurements weight
78+
public static double[][][][] MEASWEIGHTS = null;
7779
public static double[] MEASWEIGHT = { 1, 1, 1, 1, 1, 1, 1,
7880
1, 1, 1, 1, 1, 1,
7981
1, 1, 1, 1, 1, 1,

‎dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Histo.java‎

Lines changed: 126 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -516,9 +516,8 @@ private void processEvent(Event event, Event shifted) {
516516
else
517517
this.leftright[sector-1][it][ip].getH1F("hi-lL" + hit.layer).fill(hit.time);
518518
}
519-
}
520-
this.vertex[it][ip].getH1F("hi-S" + sector).fill(vz);
521-
519+
}
520+
this.vertex[it][ip].getH1F("hi-S" + sector).fill(vz);
522521
}
523522
}
524523
}
@@ -650,6 +649,9 @@ public void analyzeHisto(int fit, int vertexFit) {
650649
// if(l>24) this.parErrors[is][it][ip][l] *= 2;
651650
}
652651
}
652+
else {
653+
Constants.MEASWEIGHTS[is][it][ip][l]=0;
654+
}
653655
System.out.print("\r");
654656
}
655657
H1F hvtx = vertex[it][ip].getH1F("hi-S"+s);
@@ -669,18 +671,29 @@ public void analyzeHisto(int fit, int vertexFit) {
669671
if(hvtx.getFunction().parameter(i).name().equals("sc")) isc = i;
670672
if(hvtx.getFunction().parameter(i).name().equals("scw")) iscw = i;
671673
}
672-
if(isc>=0 && iscw>=0 && hvtx.getFunction().getParameter(isc)>10) {
674+
if(isc>=0 && iscw>=0 && hvtx.getFunction().getParameter(isc)>0) {
673675
this.parValues[is][it][ip][nLayer+nTarget-1] = (hvtx.getFunction().getParameter(iscw)-Constants.SCEXIT)*Constants.SCALE;
674676
this.parErrors[is][it][ip][nLayer+nTarget-1] = Math.max(hvtx.getFunction().parameter(iscw).error()*Constants.SCALE, Constants.SCALE*dx);
675677
this.parSigmas[is][it][ip][nLayer+nTarget-1] = hvtx.getFunction().getParameter(2)*Constants.SCALE;
676678
}
677-
if(itl>=0 && hvtx.getFunction().getParameter(0)>10) {
679+
else {
680+
Constants.MEASWEIGHTS[is][it][ip][nLayer+nTarget-1]=0;
681+
}
682+
if(itl>=0 && hvtx.getFunction().getParameter(0)>0) {
678683
this.parValues[is][it][ip][nLayer+nTarget-2] = (hvtx.getFunction().getParameter(itl)-Constants.TARGETLENGTH)*Constants.SCALE;
679684
this.parErrors[is][it][ip][nLayer+nTarget-2] = Math.max(hvtx.getFunction().parameter(itl).error()*Constants.SCALE, Constants.SCALE*dx);
680685
this.parSigmas[is][it][ip][nLayer+nTarget-2] = hvtx.getFunction().getParameter(2)*Constants.SCALE;
681686
}
687+
else {
688+
Constants.MEASWEIGHTS[is][it][ip][nLayer+nTarget-2]=0;
689+
}
682690
}
683691
}
692+
else {
693+
Constants.MEASWEIGHTS[is][it][ip][0]=0;
694+
Constants.MEASWEIGHTS[is][it][ip][nLayer+nTarget-2]=0;
695+
Constants.MEASWEIGHTS[is][it][ip][nLayer+nTarget-1]=0;
696+
}
684697
}
685698
}
686699
}
@@ -714,12 +727,14 @@ public void getFailedFitStats() {
714727
!residuals[is][it][ip].getH1F("hi-L"+l).getFunction().isFitValid()) {
715728
nfailed++;
716729
LOGGER.log(Level.WARNING, String.format("\tResidual fit for sector=%1d theta bin=%1d phi bin=%1d layer=%2d FAILED",is+1,it,ip,l));
730+
residuals[is][it][ip].getH1F("hi-L"+l).getFunction().setChiSquare(-Math.abs(residuals[is][it][ip].getH1F("hi-L"+l).getFunction().getChiSquare()));
717731
}
718732
}
719733
if(vertex[it][ip].getH1F("hi-S"+(is+1)).getFunction()!=null &&
720734
!vertex[it][ip].getH1F("hi-S"+(is+1)).getFunction().isFitValid()) {
721735
nfailed++;
722736
LOGGER.log(Level.WARNING, String.format("\tVertex fit for sector=%1d theta bin=%1d phi bin=%1d FAILED",is+1,it,ip));
737+
vertex[it][ip].getH1F("hi-S"+(is+1)).getFunction().setChiSquare(-Math.abs(vertex[it][ip].getH1F("hi-S"+(is+1)).getFunction().getChiSquare()));
723738
}
724739
}
725740
}
@@ -756,7 +771,7 @@ public double[] getParValues(String parameter, int sector, int itheta, int iphi)
756771
else if(parameter.equals("LR"))
757772
return this.lrValues[sector-1][itheta][iphi];
758773
else
759-
return this.parValues[sector-1][itheta][iphi];
774+
return this.getWeightedParValues(sector, itheta, iphi);
760775
}
761776

762777
public double[] getParErrors(String parameter, int sector, int itheta, int iphi) {
@@ -771,7 +786,7 @@ public double[] getParErrors(String parameter, int sector, int itheta, int iphi)
771786
else if(parameter.equals("LR"))
772787
return this.zeroes[sector-1][itheta][iphi];
773788
else
774-
return this.parErrors[sector-1][itheta][iphi];
789+
return this.getWeightedParErrors(sector, itheta, iphi);
775790
}
776791

777792
public double[] getParSigmas(String parameter, int sector, int itheta, int iphi) {
@@ -789,6 +804,22 @@ else if(parameter.equals("LR"))
789804
return this.parSigmas[sector-1][itheta][iphi];
790805
}
791806

807+
private double[] getWeightedParValues(int sector, int itheta, int iphi) {
808+
double[] wpars = new double[parValues[sector-1][itheta][iphi].length];
809+
for(int il=0; il<wpars.length; il++) {
810+
wpars[il] = parValues[sector-1][itheta][iphi][il]*Constants.MEASWEIGHTS[sector-1][itheta][iphi][il];
811+
}
812+
return wpars;
813+
}
814+
815+
private double[] getWeightedParErrors(int sector, int itheta, int iphi) {
816+
double[] wpars = new double[parErrors[sector-1][itheta][iphi].length];
817+
for(int il=0; il<wpars.length; il++) {
818+
wpars[il] = parErrors[sector-1][itheta][iphi][il]*Constants.MEASWEIGHTS[sector-1][itheta][iphi][il];
819+
}
820+
return wpars;
821+
}
822+
792823
public EmbeddedCanvasTabbed getElectronPlots() {
793824
EmbeddedCanvasTabbed canvas = new EmbeddedCanvasTabbed("electron", "binning", "offset");
794825
canvas.getCanvas("electron").draw(electron);
@@ -903,14 +934,27 @@ private static boolean fitVertex(int mode, H1F histo) {
903934
case 8:
904935
Histo.fitRGAS18Vertex(histo);
905936
break;
937+
case 9:
938+
Histo.fitRGKVertex(histo);
939+
break;
906940
default:
907-
if(histo.getFunction()!=null)
908-
histo.getFunction().setFitValid(true);
941+
if(histo.getFunction()!=null) {
942+
if(histo.getFunction().getChiSquare()>0)
943+
histo.getFunction().setFitValid(true);
944+
else
945+
histo.getFunction().setFitValid(false);
946+
}
909947
break;
910948
}
911-
histo.getFunction().setStatBoxFormat("%.2f");
912-
histo.getFunction().setStatBoxErrorFormat("%.2f");
913-
return histo.getFunction().isFitValid();
949+
if(histo.getFunction()!=null) {
950+
if(histo.getFunction().getChiSquare()/histo.getFunction().getNDF()>Constants.CHI2MAX)
951+
histo.getFunction().setFitValid(false);
952+
histo.getFunction().setStatBoxFormat("%.2f");
953+
histo.getFunction().setStatBoxErrorFormat("%.2f");
954+
return histo.getFunction().isFitValid();
955+
}
956+
else
957+
return false;
914958
}
915959

916960

@@ -992,7 +1036,10 @@ else if(fit>0) {
9921036
}
9931037
}
9941038
else if(histo.getFunction()!=null) {
995-
histo.getFunction().setFitValid(true);
1039+
if(histo.getFunction().getChiSquare()>0)
1040+
histo.getFunction().setFitValid(true);
1041+
else
1042+
histo.getFunction().setFitValid(false);
9961043
}
9971044
histo.getFunction().setStatBoxFormat("%.1f");
9981045
histo.getFunction().setStatBoxErrorFormat("%.1f");
@@ -1384,6 +1431,72 @@ public static void fitRGDVertex(H1F histo) {
13841431
}
13851432

13861433
/**
1434+
<<<<<<< HEAD
1435+
* 3-peaks vertex fitting function
1436+
* Peaks correspond to: target windows and scattering chamber exit window
1437+
* Initialized according to:
1438+
* - chosen target length (TARGETLENGTH),
1439+
* - target exit window position (TARGETPOS)
1440+
* - distance between target exit window and insulation foil (WINDOWDIST)
1441+
* - distance between the scattering chamber exit window and the target center (SCEXIT)
1442+
* Includes two wide Gaussians to account for target residual gas
1443+
* and the air outside the scattering chamber
1444+
* @param histo
1445+
*/
1446+
public static void fitRGKVertex(H1F histo) {
1447+
int nbin = histo.getData().length;
1448+
double dx = histo.getDataX(1)-histo.getDataX(0);
1449+
//find downstream window
1450+
int ibin0 = Histo.getMaximumBinBetween(histo, histo.getDataX(0), (Constants.TARGETPOS+Constants.SCEXIT)/2);
1451+
//check if the found maximum is the first or second peak, ibin is tentative upstream window
1452+
int ibin1 = Math.max(0, ibin0 - (int)(Constants.TARGETLENGTH/dx));
1453+
int ibin2 = Math.min(nbin-1, ibin0 + (int)(Constants.TARGETLENGTH/dx));
1454+
if(histo.getBinContent(ibin1)<histo.getBinContent(ibin2)) {
1455+
ibin1 = ibin0;
1456+
ibin0 = ibin2;
1457+
}
1458+
int ibinsc = Histo.getMaximumBinBetween(histo, (Constants.SCEXIT+Constants.TARGETCENTER)*0.8, (Constants.SCEXIT+Constants.TARGETCENTER)*1.2);
1459+
1460+
double mean = histo.getDataX(ibin0);
1461+
double amp = histo.getBinContent(ibin0);
1462+
double sigma = 0.5;
1463+
double sc = histo.getBinContent(ibinsc);
1464+
double air = histo.getBinContent(ibinsc + ((int) (sigma*6/dx)));
1465+
double scw = Constants.SCEXIT;
1466+
if(sc>10) scw = histo.getDataX(ibinsc)-mean+Constants.TARGETLENGTH/2;
1467+
double bg = histo.getBinContent((ibin1+ibin0)/2);
1468+
String function = "[ampU]*gaus(x,[exw]-[tl],[sigma])+"
1469+
+ "[ampD]*gaus(x,[exw],[sigma])+"
1470+
+ "[bg]*gaus(x,[exw]-[tl]/2,[tl]*0.6)+"
1471+
+ "[sc]*gaus(x,[exw]+[scw]-[tl]/2,[sigma])+"
1472+
+ "[air]*gaus(x,[exw]+[scw]-[tl]/2+[adelta],[asigma])";
1473+
F1D f1_vtx = new F1D("f"+histo.getName(), function, -10, 10);
1474+
f1_vtx.setLineColor(2);
1475+
f1_vtx.setLineWidth(2);
1476+
f1_vtx.setOptStat("1111111111111");
1477+
f1_vtx.setParameter(0, amp);
1478+
f1_vtx.setParameter(1, mean);
1479+
f1_vtx.setParameter(2, Constants.TARGETLENGTH);
1480+
f1_vtx.setParLimits(2, Constants.TARGETLENGTH*0.9, Constants.TARGETLENGTH*1.1);
1481+
f1_vtx.setParameter(3, sigma);
1482+
f1_vtx.setParameter(4, amp);
1483+
f1_vtx.setParameter(5, bg);
1484+
f1_vtx.setParameter(6, sc);
1485+
f1_vtx.setParameter(7, scw);
1486+
f1_vtx.setParLimits(7, (Constants.SCEXIT)*0.7, (Constants.SCEXIT)*1.3);
1487+
f1_vtx.setParameter(8, air);
1488+
f1_vtx.setParameter(9, sigma*3);
1489+
f1_vtx.setParLimits(9, 0, sigma*8);
1490+
f1_vtx.setParameter(10, sigma*8);
1491+
f1_vtx.setRange(mean-Constants.TARGETLENGTH*2.0,Constants.SCEXIT+Constants.TARGETLENGTH*0.6);
1492+
// histo.setFunction(f1_vtx);
1493+
DataFitter.fit(f1_vtx, histo, "Q"); //No options uses error for sigma
1494+
// if(f1_vtx.getParameter(6)<f1_vtx.getParameter(0)/4) f1_vtx.setParameter(6, 0);
1495+
}
1496+
1497+
/**
1498+
=======
1499+
>>>>>>> master
13871500
* 4-peaks vertex fitting function
13881501
* Peaks correspond to: target windows and scattering chamber exit window
13891502
* Initialized according to:

‎dc/utilities/createSkims.csh‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ else
2424
set vars = "r0 r1_x r1_y r1_z r1_cy r1_cz r2_x r2_y r2_z r2_cy r2_cz r3_x r3_y r3_z r3_cy r3_cz"
2525
endif
2626

27-
set schema = `ls $indir/r0`
27+
set schema = `ls $indir/r0/*/recon/README.json | awk -F"$indir/r0/" '{print $2}' | awk -F"/" '{print $1}'`
2828

2929
echo
3030
echo reading input files from $indir with subdirectory $schema

‎dc/utilities/kinematics.groovy‎

Lines changed: 28 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -46,11 +46,24 @@ public class Kinematics {
4646
private Map<String, DataGroup> histos = new LinkedHashMap<>();
4747

4848

49-
public Kinematics(double energy) {
49+
public Kinematics(double energy, String vertex) {
5050
this.initGraphics();
51+
this.initVertexPar(vertex);
5152
this.setEbeam(energy);
5253
}
5354

55+
private void initVertexPar(String pars) {
56+
if(!pars.isEmpty()) {
57+
double[] parValues = new double[pars.split(":").length];
58+
for(int i=0; i<parValues.length; i++) {
59+
parValues[i] = Double.parseDouble(pars.split(":")[i]);
60+
}
61+
if(parValues.length>0) targetPos = parValues[0];
62+
if(parValues.length>1) targetLength = parValues[1];
63+
if(parValues.length>2) scWindow = parValues[2]+targetPos;
64+
}
65+
}
66+
5467
public double getEbeam() {
5568
return ebeam;
5669
}
@@ -724,12 +737,22 @@ public class Kinematics {
724737
parser.getOptionParser("-process").addOption("-beam" ,"10.6", "beam energy in GeV");
725738
parser.getOptionParser("-process").addOption("-display" ,"1", "display histograms (0/1)");
726739
parser.getOptionParser("-process").addOption("-stats" ,"", "histogram stat option");
740+
parser.getOptionParser("-process").addOption("-vertpar" , "", "comma-separated vertex function parameters, default values are for Spring19 cryotarget with:\n" +
741+
"\t\t- -3.5: target cell exit window position,\n" +
742+
"\t\t- 5.0: target length,\n" +
743+
"\t\t- 27.3: distance between the scattering chamber exit window and the target center,\n" +
744+
"\t\t leave empty to use defaults; units are cm");
727745

728746
// valid options for histogram-base analysis
729747
parser.addCommand("-plot", "plot histogram files");
730748
parser.getOptionParser("-plot").addOption("-beam" ,"10.6", "beam energy in GeV");
731749
parser.getOptionParser("-plot").addOption("-display" ,"1", "display histograms (0/1)");
732750
parser.getOptionParser("-plot").addOption("-stats" ,"", "set histogram stat option");
751+
parser.getOptionParser("-plot").addOption("-vertpar" , "", "comma-separated vertex function parameters, default values are for Spring19 cryotarget with:\n" +
752+
"\t\t- -3.5: target cell exit window position,\n" +
753+
"\t\t- 5.0: target length,\n" +
754+
"\t\t- 27.3: distance between the scattering chamber exit window and the target center,\n" +
755+
"\t\t leave empty to use defaults; units are cm");
733756

734757
parser.parse(args);
735758

@@ -742,6 +765,7 @@ public class Kinematics {
742765
if(parser.getCommand().equals("-process")) {
743766
int maxEvents = parser.getOptionParser("-process").getOption("-nevent").intValue();
744767
double beamEnergy = parser.getOptionParser("-process").getOption("-beam").doubleValue();
768+
String vertexPar = parser.getOptionParser("-process").getOption("-vertpar").stringValue();
745769
String namePrefix = parser.getOptionParser("-process").getOption("-o").stringValue();
746770
String histoName = "histo.hipo";
747771
if(!namePrefix.isEmpty()) {
@@ -758,7 +782,7 @@ public class Kinematics {
758782
System.exit(0);
759783
}
760784

761-
analysis = new Kinematics(beamEnergy);
785+
analysis = new Kinematics(beamEnergy, vertexPar);
762786

763787
ProgressPrintout progress = new ProgressPrintout();
764788

@@ -790,6 +814,7 @@ public class Kinematics {
790814

791815
if(parser.getCommand().equals("-plot")) {
792816
double beamEnergy = parser.getOptionParser("-plot").getOption("-beam").doubleValue();
817+
String vertexPar = parser.getOptionParser("-plot").getOption("-vertpar").stringValue();
793818
optStats = parser.getOptionParser("-plot").getOption("-stats").stringValue();
794819
openWindow = parser.getOptionParser("-plot").getOption("-display").intValue()!=0;
795820
if(!openWindow) System.setProperty("java.awt.headless", "true");
@@ -801,7 +826,7 @@ public class Kinematics {
801826
System.exit(0);
802827
}
803828

804-
analysis = new Kinematics(beamEnergy);
829+
analysis = new Kinematics(beamEnergy, vertexPar);
805830
analysis.readHistos(inputList.get(0));
806831
analysis.analyzeHistos();
807832
}

0 commit comments

Comments
 (0)