14 #include "TDirectory.h"
19 #include <JANA/JApplication.h>
20 #include <JANA/JFactory.h>
51 milleWriter =
new Mille(
"fieldon_mille_out.mil");
53 gDirectory->mkdir(
"AlignmentConstants");
54 gDirectory->cd(
"AlignmentConstants");
57 HistCurrentConstantsCDC =
new TProfile(
"CDCAlignmentConstants",
"Constants Used for CDC Alignment (In MILLEPEDE Order)", 20000 ,0.5, 20000.5);
58 HistCurrentConstantsFDC =
new TProfile(
"FDCAlignmentConstants",
"Constants Used for FDC Alignment (In MILLEPEDE Order)", 26000 ,0.5, 26000.5);
78 jerr <<
" ********No Field******** Plugin MilleFieldOn is for use with magnetic field!!! Aborting " << 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(
"CDC/timing_offsets",vals)==
false){
108 for(
unsigned int i=0; i<vals.size(); i++){
109 map<string,double> &row = vals[i];
111 HistCurrentConstantsCDC->Fill(16000+(i+1),row[
"t0"]);
115 if (jcalib->Get(
"FDC/cathode_alignment",vals)==
false){
116 for(
unsigned int i=0; i<vals.size(); i++){
117 map<string,double> &row = vals[i];
119 HistCurrentConstantsFDC->Fill((i+1)*1000+101,row[
"dU"]);
120 HistCurrentConstantsFDC->Fill((i+1)*1000+102,row[
"dV"]);
121 HistCurrentConstantsFDC->Fill((i+1)*1000+103,row[
"dPhiU"]);
122 HistCurrentConstantsFDC->Fill((i+1)*1000+104,row[
"dPhiV"]);
126 if (jcalib->Get(
"FDC/cell_offsets",vals)==
false){
127 for(
unsigned int i=0; i<vals.size(); i++){
128 map<string,double> &row = vals[i];
130 HistCurrentConstantsFDC->Fill((i+1)*1000+1,row[
"xshift"]);
131 HistCurrentConstantsFDC->Fill((i+1)*1000+100,row[
"yshift"]);
135 if (jcalib->Get(
"FDC/cell_rotations",vals)==
false){
136 for(
unsigned int i=0; i<vals.size(); i++){
137 map<string,double> &row = vals[i];
138 HistCurrentConstantsFDC->Fill((i+1)*1000+2,row[
"dPhiX"]);
139 HistCurrentConstantsFDC->Fill((i+1)*1000+3,row[
"dPhiY"]);
140 HistCurrentConstantsFDC->Fill((i+1)*1000+4,row[
"dPhiZ"]);
143 if (jcalib->Get(
"FDC/wire_alignment",vals)==
false){
144 for(
unsigned int i=0; i<vals.size(); i++){
145 map<string,double> &row = vals[i];
147 HistCurrentConstantsFDC->Fill((i+1)*1000+5,row[
"dZ"]);
151 if (jcalib->Get(
"FDC/package1/strip_gains_v2",vals)==
false){
152 for(
unsigned int i=0; i<vals.size(); i++){
153 map<string,double> &row = vals[i];
156 for (
unsigned int j=1; j<= 216; j++)
160 HistCurrentConstantsFDC->Fill((i/2+1)*1000+600+j,row[name]);
163 HistCurrentConstantsFDC->Fill((i/2+1)*1000+300+j,row[name]);
169 if (jcalib->Get(
"FDC/package2/strip_gains_v2",vals)==
false){
170 for(
unsigned int i=0; i<vals.size(); i++){
171 map<string,double> &row = vals[i];
174 for (
unsigned int j=1; j<= 216; j++)
178 HistCurrentConstantsFDC->Fill((i/2+7)*1000+600+j,row[name]);
181 HistCurrentConstantsFDC->Fill((i/2+7)*1000+300+j,row[name]);
186 if (jcalib->Get(
"FDC/package3/strip_gains_v2",vals)==
false){
187 for(
unsigned int i=0; i<vals.size(); i++){
188 map<string,double> &row = vals[i];
191 for (
unsigned int j=1; j<= 216; j++)
195 HistCurrentConstantsFDC->Fill((i/2+13)*1000+600+j,row[name]);
198 HistCurrentConstantsFDC->Fill((i/2+13)*1000+300+j,row[name]);
203 if (jcalib->Get(
"FDC/package4/strip_gains_v2",vals)==
false){
204 for(
unsigned int i=0; i<vals.size(); i++){
205 map<string,double> &row = vals[i];
208 for (
unsigned int j=1; j<= 216; j++)
212 HistCurrentConstantsFDC->Fill((i/2+19)*1000+600+j,row[name]);
215 HistCurrentConstantsFDC->Fill((i/2+19)*1000+300+j,row[name]);
220 if (jcalib->Get(
"FDC/package1/wire_timing_offsets",vals)==
false){
221 for(
unsigned int i=0; i<vals.size(); i++){
222 map<string,double> &row = vals[i];
225 for (
unsigned int j=1; j<= 96; j++)
228 HistCurrentConstantsFDC->Fill((i+1)*1000+900+j,row[name]);
232 if (jcalib->Get(
"FDC/package2/wire_timing_offsets",vals)==
false){
233 for(
unsigned int i=0; i<vals.size(); i++){
234 map<string,double> &row = vals[i];
237 for (
unsigned int j=1; j<= 96; j++)
240 HistCurrentConstantsFDC->Fill((i+7)*1000+900+j,row[name]);
244 if (jcalib->Get(
"FDC/package3/wire_timing_offsets",vals)==
false){
245 for(
unsigned int i=0; i<vals.size(); i++){
246 map<string,double> &row = vals[i];
249 for (
unsigned int j=1; j<= 96; j++)
252 HistCurrentConstantsFDC->Fill((i+13)*1000+900+j,row[name]);
256 if (jcalib->Get(
"FDC/package4/wire_timing_offsets",vals)==
false){
257 for(
unsigned int i=0; i<vals.size(); i++){
258 map<string,double> &row = vals[i];
261 for (
unsigned int j=1; j<= 96; j++)
264 HistCurrentConstantsFDC->Fill((i+19)*1000+900+j,row[name]);
270 if (jcalib->Get(
"FDC/strip_pitches_v2",vals)==
false){
271 for(
unsigned int i=0; i<vals.size(); i++){
272 map<string,double> &row = vals[i];
273 HistCurrentConstantsFDC->Fill((i+1)*1000 + 200, row[
"U_SP_1"]);
274 HistCurrentConstantsFDC->Fill((i+1)*1000 + 201, row[
"U_G_1"]);
275 HistCurrentConstantsFDC->Fill((i+1)*1000 + 202, row[
"U_SP_2"]);
276 HistCurrentConstantsFDC->Fill((i+1)*1000 + 203, row[
"U_G_2"]);
277 HistCurrentConstantsFDC->Fill((i+1)*1000 + 204, row[
"U_SP_3"]);
278 HistCurrentConstantsFDC->Fill((i+1)*1000 + 205, row[
"V_SP_1"]);
279 HistCurrentConstantsFDC->Fill((i+1)*1000 + 206, row[
"V_G_1"]);
280 HistCurrentConstantsFDC->Fill((i+1)*1000 + 207, row[
"V_SP_2"]);
281 HistCurrentConstantsFDC->Fill((i+1)*1000 + 208, row[
"V_G_2"]);
282 HistCurrentConstantsFDC->Fill((i+1)*1000 + 209, row[
"V_SP_3"]);
294 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};
296 vector<const DChargedTrack *> trackVector;
297 loop->Get(trackVector);
299 for (
size_t i = 0; i < trackVector.size(); i++){
303 if (bestHypothesis == NULL)
continue;
306 double trackingFOM = TMath::Prob(thisTimeBasedTrack->chisq, thisTimeBasedTrack->Ndof);
310 float trackingFOMCut = 0.001;
311 if(trackingFOM < trackingFOMCut)
continue;
314 if (!thisTimeBasedTrack->IsSmoothed)
continue;
316 vector<DTrackFitter::pull_t> pulls = thisTimeBasedTrack->pulls;
322 for (
size_t iPull = 0; iPull < pulls.size(); iPull++){
325 float err = pulls[iPull].err;
326 float errc = pulls[iPull].errc;
327 if (cdc_hit == NULL) {pullsNDF+=2;nFDC++;isCDCOnly=
false;}
328 else {pullsNDF++;nCDC++;}
329 if ( err != err || errc != errc) anyNaN=
true;
333 if (anyNaN)
continue;
334 if(pullsNDF != thisTimeBasedTrack->Ndof)
continue;
336 int trackingNDFCut = 22;
337 if (isCDCOnly) trackingNDFCut = 10;
339 if( thisTimeBasedTrack->Ndof < trackingNDFCut)
continue;
348 japp->RootWriteLock();
349 for (
size_t iPull = 0; iPull < pulls.size(); iPull++){
352 float resi = pulls[iPull].resi;
353 float err = pulls[iPull].err;
358 const DFDCPseudo *fdc_hit = pulls[iPull].fdc_hit;
362 float resic = pulls[iPull].resic;
363 float errc = pulls[iPull].errc;
365 vector<double> trackDerivatives = pulls[iPull].trackDerivatives;
366 if (trackDerivatives.size()==0)
continue;
368 if (fdc_hit != NULL && fdc_hit->
status == 6) {
376 unsigned int layerOffset = 100000 + thisHit->
wire->
layer * 1000;
381 float localDerW[NLC];
382 float globalDerW[NGL_W];
395 milleWriter->mille(NLC, localDerW, NGL_W, globalDerW, labelW, resi, err);
398 const int NGL_C = 26;
399 float localDerC[NLC];
400 float globalDerC[NGL_C];
406 globalDerC[0] = -1.0; labelC[0] = layerOffset + 100;
420 globalDerC[10] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[0]; labelC[10] = layerOffset + 200;
421 globalDerC[11] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[1]; labelC[11] = layerOffset + 201;
422 globalDerC[12] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[2]; labelC[12] = layerOffset + 202;
423 globalDerC[13] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[3]; labelC[13] = layerOffset + 203;
424 globalDerC[14] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripPitchDerivatives[4]; labelC[14] = layerOffset + 204;
425 globalDerC[15] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[5]; labelC[15] = layerOffset + 205;
426 globalDerC[16] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[6]; labelC[16] = layerOffset + 206;
427 globalDerC[17] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[7]; labelC[17] = layerOffset + 207;
428 globalDerC[18] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[8]; labelC[18] = layerOffset + 208;
429 globalDerC[19] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripPitchDerivatives[9]; labelC[19] = layerOffset + 209;
432 vector<const DFDCCathodeCluster *> cathodeClusters;
433 thisHit->Get(cathodeClusters);
434 unsigned int gainLabels[6]={};
435 for(
unsigned int j = 0; j< cathodeClusters.size(); j++){
436 if(cathodeClusters[j]->plane == 3) {
438 for (; k < cathodeClusters[j]->members.size(); k++){
439 if (thisHit->
cluster_u.
index(0) == cathodeClusters[j]->members[k]->element)
break;
445 else if (cathodeClusters[j]->plane == 1) {
447 for (; k < cathodeClusters[j]->members.size(); k++){
448 if (thisHit->
cluster_v.
index(0) == cathodeClusters[j]->members[k]->element)
break;
456 globalDerC[20] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripGainDerivatives[0]; labelC[20] = layerOffset + 300 + gainLabels[0];
457 globalDerC[21] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripGainDerivatives[1]; labelC[21] = layerOffset + 300 + gainLabels[1];
458 globalDerC[22] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaU]*fdcStripGainDerivatives[2]; labelC[22] = layerOffset + 300 + gainLabels[2];
459 globalDerC[23] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripGainDerivatives[3]; labelC[23] = layerOffset + 600 + gainLabels[3];
460 globalDerC[24] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripGainDerivatives[4]; labelC[24] = layerOffset + 600 + gainLabels[4];
461 globalDerC[25] = -pseudoAlignmentDerivatives[
FDCPseudoD::dSddeltaV]*fdcStripGainDerivatives[5]; labelC[25] = layerOffset + 600 + gainLabels[5];
463 milleWriter->mille(NLC, localDerC, NGL_C, globalDerC, labelC, resic, errc);
466 if (cdc_hit != NULL){
477 float globalDer[NGL];
485 globalDer[0]=0.0; label[0]=1;
486 globalDer[1]=0.0; label[0]=2;
487 globalDer[2]=0.0; label[0]=3;
488 globalDer[3]=0.0; label[0]=4;
489 globalDer[4]=0.0; label[0]=5;
490 globalDer[5]=0.0; label[0]=6;
547 label[6]=1000 + (straw_offset[thisWire->
ring]+(thisWire->
straw-1))*4 + 1;
550 jout <<
" Dumping deltaXu derivatives============ Wire " << thisWire->
ring <<
" Straw " << thisWire->
straw << endl;
551 jout <<
" Total = " << globalDer[6] << endl;
566 label[7]=1000 + (straw_offset[thisWire->
ring]+(thisWire->
straw-1))*4 + 2;
574 label[8]=1000 + (straw_offset[thisWire->
ring]+(thisWire->
straw-1))*4 + 3;
582 label[9]=1000 + (straw_offset[thisWire->
ring]+(thisWire->
straw-1))*4 + 4;
585 label[10]=16000+straw_offset[thisWire->
ring]+thisWire->
straw;
587 milleWriter->mille(NLC, localDer, NGL, globalDer, label, resi, err);
int status
status word for pseudopoint
vector< double > GetCDCWireDerivatives()
sprintf(text,"Post KinFit Cut")
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
const DTrackTimeBased * Get_TrackTimeBased(void) const
const DFDCWire * wire
DFDCWire for this wire.
class DFDCPseudo: definition for a reconstructed point in the FDC
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
vector< double > GetFDCPseudoAlignmentDerivatives()
centroid_t cluster_v
Cathode cluster centroids, Useful for gain/strip pitch calibration.
vector< double > GetFDCStripGainDerivatives()
JEventProcessor_MilleFieldOn()
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)
const vector< double > GetFDCStripPitchDerivatives()
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
~JEventProcessor_MilleFieldOn()
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
jerror_t init(void)
Called once at program start.
jerror_t fini(void)
Called after last event of last event source has been processed.