13 #include "TDirectory.h"
18 #include <JANA/JApplication.h>
19 #include <JANA/JFactory.h>
50 milleWriter =
new Mille(
"nofield_mille_out.mil");
52 gDirectory->mkdir(
"AlignmentConstants");
53 gDirectory->cd(
"AlignmentConstants");
56 HistCurrentConstantsCDC =
new TProfile(
"CDCAlignmentConstants",
"Constants Used for CDC Alignment (In MILLEPEDE Order)", 16000 ,0.5, 16000.5);
57 HistCurrentConstantsFDC =
new TProfile(
"FDCAlignmentConstants",
"Constants Used for FDC Alignment (In MILLEPEDE Order)", 26000 ,0.5, 26000.5);
77 jerr <<
" Plugin FDC_MilleFieldOff Must be run with zero magnetic field!!! Aborting " << endl;
78 jerr <<
" Use -PBFIELD_TYPE=NoField " << endl;
83 JCalibration * jcalib = eventLoop->GetJCalibration();
84 vector<map<string,double> >vals;
85 if (jcalib->Get(
"CDC/global_alignment",vals)==
false){
86 map<string,double> &row = vals[0];
88 HistCurrentConstantsCDC->Fill(1,row[
"dX"]);
89 HistCurrentConstantsCDC->Fill(2,row[
"dY"]);
90 HistCurrentConstantsCDC->Fill(3,row[
"dZ"]);
91 HistCurrentConstantsCDC->Fill(4,row[
"dPhiX"]);
92 HistCurrentConstantsCDC->Fill(5,row[
"dPhiY"]);
93 HistCurrentConstantsCDC->Fill(6,row[
"dPhiZ"]);
96 if (jcalib->Get(
"CDC/wire_alignment",vals)==
false){
97 for(
unsigned int i=0; i<vals.size(); i++){
98 map<string,double> &row = vals[i];
100 HistCurrentConstantsCDC->Fill(1000+(i*4+1),row[
"dxu"]);
101 HistCurrentConstantsCDC->Fill(1000+(i*4+2),row[
"dyu"]);
102 HistCurrentConstantsCDC->Fill(1000+(i*4+3),row[
"dxd"]);
103 HistCurrentConstantsCDC->Fill(1000+(i*4+4),row[
"dyd"]);
107 if (jcalib->Get(
"FDC/cathode_alignment",vals)==
false){
108 for(
unsigned int i=0; i<vals.size(); i++){
109 map<string,double> &row = vals[i];
111 HistCurrentConstantsFDC->Fill((i+1)*1000+101,row[
"dU"]);
112 HistCurrentConstantsFDC->Fill((i+1)*1000+102,row[
"dV"]);
113 HistCurrentConstantsFDC->Fill((i+1)*1000+103,row[
"dPhiU"]);
114 HistCurrentConstantsFDC->Fill((i+1)*1000+104,row[
"dPhiV"]);
118 if (jcalib->Get(
"FDC/cell_offsets",vals)==
false){
119 for(
unsigned int i=0; i<vals.size(); i++){
120 map<string,double> &row = vals[i];
122 HistCurrentConstantsFDC->Fill((i+1)*1000+1,row[
"xshift"]);
123 HistCurrentConstantsFDC->Fill((i+1)*1000+100,row[
"yshift"]);
127 if (jcalib->Get(
"FDC/cell_rotations",vals)==
false){
128 for(
unsigned int i=0; i<vals.size(); i++){
129 map<string,double> &row = vals[i];
130 HistCurrentConstantsFDC->Fill((i+1)*1000+2,row[
"dPhiX"]);
131 HistCurrentConstantsFDC->Fill((i+1)*1000+3,row[
"dPhiY"]);
132 HistCurrentConstantsFDC->Fill((i+1)*1000+4,row[
"dPhiZ"]);
135 if (jcalib->Get(
"FDC/wire_alignment",vals)==
false){
136 for(
unsigned int i=0; i<vals.size(); i++){
137 map<string,double> &row = vals[i];
139 HistCurrentConstantsFDC->Fill((i+1)*1000+5,row[
"dZ"]);
143 if (jcalib->Get(
"FDC/package1/strip_gains_v2",vals)==
false){
144 for(
unsigned int i=0; i<vals.size(); i++){
145 map<string,double> &row = vals[i];
148 for (
unsigned int j=1; j<= 216; j++)
152 HistCurrentConstantsFDC->Fill((i/2+1)*1000+600+j,row[name]);
155 HistCurrentConstantsFDC->Fill((i/2+1)*1000+300+j,row[name]);
161 if (jcalib->Get(
"FDC/package2/strip_gains_v2",vals)==
false){
162 for(
unsigned int i=0; i<vals.size(); i++){
163 map<string,double> &row = vals[i];
166 for (
unsigned int j=1; j<= 216; j++)
170 HistCurrentConstantsFDC->Fill((i/2+7)*1000+600+j,row[name]);
173 HistCurrentConstantsFDC->Fill((i/2+7)*1000+300+j,row[name]);
178 if (jcalib->Get(
"FDC/package3/strip_gains_v2",vals)==
false){
179 for(
unsigned int i=0; i<vals.size(); i++){
180 map<string,double> &row = vals[i];
183 for (
unsigned int j=1; j<= 216; j++)
187 HistCurrentConstantsFDC->Fill((i/2+13)*1000+600+j,row[name]);
190 HistCurrentConstantsFDC->Fill((i/2+13)*1000+300+j,row[name]);
195 if (jcalib->Get(
"FDC/package4/strip_gains_v2",vals)==
false){
196 for(
unsigned int i=0; i<vals.size(); i++){
197 map<string,double> &row = vals[i];
200 for (
unsigned int j=1; j<= 216; j++)
204 HistCurrentConstantsFDC->Fill((i/2+19)*1000+600+j,row[name]);
207 HistCurrentConstantsFDC->Fill((i/2+19)*1000+300+j,row[name]);
212 if (jcalib->Get(
"FDC/strip_pitches_v2",vals)==
false){
213 for(
unsigned int i=0; i<vals.size(); i++){
214 map<string,double> &row = vals[i];
215 HistCurrentConstantsFDC->Fill((i+1)*1000 + 200, row[
"U_SP_1"]);
216 HistCurrentConstantsFDC->Fill((i+1)*1000 + 201, row[
"U_G_1"]);
217 HistCurrentConstantsFDC->Fill((i+1)*1000 + 202, row[
"U_SP_2"]);
218 HistCurrentConstantsFDC->Fill((i+1)*1000 + 203, row[
"U_G_2"]);
219 HistCurrentConstantsFDC->Fill((i+1)*1000 + 204, row[
"U_SP_3"]);
220 HistCurrentConstantsFDC->Fill((i+1)*1000 + 205, row[
"V_SP_1"]);
221 HistCurrentConstantsFDC->Fill((i+1)*1000 + 206, row[
"V_G_1"]);
222 HistCurrentConstantsFDC->Fill((i+1)*1000 + 207, row[
"V_SP_2"]);
223 HistCurrentConstantsFDC->Fill((i+1)*1000 + 208, row[
"V_G_2"]);
224 HistCurrentConstantsFDC->Fill((i+1)*1000 + 209, row[
"V_SP_3"]);
236 int straw_offset[29] = {0,0,42,84,138,192,258,324,404,484,577,670,776,882,1005,1128,1263,1398,1544,1690,1848,2006,2176,2346,2528,2710,2907,3104,3313};
238 vector<const DTrackCandidate *> trackVector;
239 loop->Get(trackVector,
"StraightLine");
241 for (
size_t i = 0; i < trackVector.size(); i++){
244 double trackingFOM = TMath::Prob(track->
chisq, track->
Ndof);
248 double trackingFOMCut = 0.001;
249 int trackingNDFCut = 18;
251 if(trackingFOM < trackingFOMCut)
continue;
252 if(track->
Ndof < trackingNDFCut)
continue;
259 vector<DTrackFitter::pull_t> pulls = track->
pulls;
260 if (pulls.size() > 75)
continue;
265 for (
size_t iPull = 0; iPull < pulls.size(); iPull++){
268 float resi = pulls[iPull].resi;
269 float err = pulls[iPull].err;
270 float resic = pulls[iPull].resic;
271 float errc = pulls[iPull].errc;
272 if (cdc_hit == NULL) isCDCOnly=
false;
273 if (resi != resi || err != err || resic != resic || errc != errc || !isfinite(resi) || !isfinite(resic)) anyNaN=
true;
275 if (fabs(resi)>10. || err > 10. || fabs(resic) > 10. || errc > 10.) anyNaN=
true;
279 if (anyNaN)
continue;
283 if(track->
Ndof < 26)
continue;
286 japp->RootWriteLock();
287 for (
size_t iPull = 0; iPull < pulls.size(); iPull++){
290 float resi = pulls[iPull].resi;
291 float err = pulls[iPull].err;
296 const DFDCPseudo *fdc_hit = pulls[iPull].fdc_hit;
300 float resic = pulls[iPull].resic;
301 float errc = pulls[iPull].errc;
303 vector<double> trackDerivatives = pulls[iPull].trackDerivatives;
304 if (trackDerivatives.size()==0) jerr <<
"Track derivative size == 0 ???" << endl;
306 if (fdc_hit != NULL && fdc_hit->
status == 6) {
314 unsigned int layerOffset = 100000 + thisHit->
wire->
layer * 1000;
319 float localDerW[NLC];
320 float globalDerW[NGL_W];
332 milleWriter->mille(NLC, localDerW, NGL_W, globalDerW, labelW, resi, err);
335 const int NGL_C = 26;
336 float localDerC[NLC];
337 float globalDerC[NGL_C];
343 globalDerC[0] = -1.0; labelC[0] = layerOffset + 100;
357 globalDerC[10] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[0]; labelC[10] = layerOffset + 200;
358 globalDerC[11] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[1]; labelC[11] = layerOffset + 201;
359 globalDerC[12] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[2]; labelC[12] = layerOffset + 202;
360 globalDerC[13] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[3]; labelC[13] = layerOffset + 203;
361 globalDerC[14] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[4]; labelC[14] = layerOffset + 204;
362 globalDerC[15] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[5]; labelC[15] = layerOffset + 205;
363 globalDerC[16] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[6]; labelC[16] = layerOffset + 206;
364 globalDerC[17] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[7]; labelC[17] = layerOffset + 207;
365 globalDerC[18] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[8]; labelC[18] = layerOffset + 208;
366 globalDerC[19] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[9]; labelC[19] = layerOffset + 209;
369 vector<const DFDCCathodeCluster *> cathodeClusters;
370 thisHit->Get(cathodeClusters);
371 unsigned int gainLabels[6]={};
372 for(
unsigned int j = 0; j< cathodeClusters.size(); j++){
373 if(cathodeClusters[j]->plane == 3) {
375 for (; k < cathodeClusters[j]->members.size(); k++){
376 if (thisHit->
cluster_u.
index(0) == cathodeClusters[j]->members[k]->element)
break;
382 else if (cathodeClusters[j]->plane == 1) {
384 for (; k < cathodeClusters[j]->members.size(); k++){
385 if (thisHit->
cluster_v.
index(0) == cathodeClusters[j]->members[k]->element)
break;
393 globalDerC[20] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripGainDerivatives[0]; labelC[20] = layerOffset + 300 + gainLabels[0];
394 globalDerC[21] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripGainDerivatives[1]; labelC[21] = layerOffset + 300 + gainLabels[1];
395 globalDerC[22] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripGainDerivatives[2]; labelC[22] = layerOffset + 300 + gainLabels[2];
396 globalDerC[23] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripGainDerivatives[3]; labelC[23] = layerOffset + 600 + gainLabels[3];
397 globalDerC[24] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripGainDerivatives[4]; labelC[24] = layerOffset + 600 + gainLabels[4];
398 globalDerC[25] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripGainDerivatives[5]; labelC[25] = layerOffset + 600 + gainLabels[5];
400 milleWriter->mille(NLC, localDerC, NGL_C, globalDerC, labelC, resic, errc);
403 if (cdc_hit != NULL){
414 float globalDer[NGL];
422 globalDer[0]=0.0; label[0]=1;
423 globalDer[1]=0.0; label[0]=2;
424 globalDer[2]=0.0; label[0]=3;
425 globalDer[3]=0.0; label[0]=4;
426 globalDer[4]=0.0; label[0]=5;
427 globalDer[5]=0.0; label[0]=6;
454 jout <<
" Dumping deltaZ derivatives============ Wire " << thisWire->
ring <<
" Straw " << thisWire->
straw << endl;
455 jout <<
" Total = " << globalDer[2] << endl;
488 jout <<
" Dumping deltaPhiZ derivatives============ Ring " << thisWire->
ring <<
" Straw " << thisWire->
straw << endl;
489 jout <<
" Total = " << globalDer[5] << endl;
504 label[6]=1000 + (straw_offset[thisWire->
ring]+(thisWire->
straw-1))*4 + 1;
507 jout <<
" Dumping deltaXu derivatives============ Wire " << thisWire->
ring <<
" Straw " << thisWire->
straw << endl;
508 jout <<
" Total = " << globalDer[6] << endl;
523 label[7]=1000 + (straw_offset[thisWire->
ring]+(thisWire->
straw-1))*4 + 2;
531 label[8]=1000 + (straw_offset[thisWire->
ring]+(thisWire->
straw-1))*4 + 3;
539 label[9]=1000 + (straw_offset[thisWire->
ring]+(thisWire->
straw-1))*4 + 4;
541 milleWriter->mille(NLC, localDer, NGL, globalDer, label, resi, err);
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
int status
status word for pseudopoint
vector< double > GetCDCWireDerivatives()
sprintf(text,"Post KinFit Cut")
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
const DFDCWire * wire
DFDCWire for this wire.
class DFDCPseudo: definition for a reconstructed point in the FDC
vector< double > GetFDCPseudoAlignmentDerivatives()
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
centroid_t cluster_v
Cathode cluster centroids, Useful for gain/strip pitch calibration.
vector< double > GetFDCStripGainDerivatives()
JEventProcessor_MilleFieldOff()
vector< DTrackFitter::pull_t > pulls
<A href="index.html#legend"> <IMG src="CORE.png" width="100"> </A>
Class to write C binary file.
void cdc_hit(Int_t &, Int_t &, Int_t &, Int_t[], Int_t, Int_t, Int_t, Int_t, Int_t, Int_t)
float chisq
Chi-squared for the track (not chisq/dof!)
const vector< double > GetFDCStripPitchDerivatives()
~JEventProcessor_MilleFieldOff()
jerror_t fini(void)
Called after last event of last event source has been processed.
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
jerror_t init(void)
Called once at program start.
int Ndof
Number of degrees of freedom in the fit.