Hall C ROOT/C++ Analyzer (hcana)
THcHallCSpectrometer.cxx
Go to the documentation of this file.
1
45 /*
46  THESE comments need to be reviewed.
47
48  th_geo is rotation about the Y axis in XZ plane to coordinates X',Y=Y',Z'
49  After th_geo rotation, ph_geo rotation is about the X' axis in the Y'Z' plane.
50
51
52  Calls THaAnalysisObject::GeoToSph to calculate the spherical angles. th_sph and ph_sph
53  th_sph is rotation about the Y axis in XZ plane to coordinates X',Y=Y',Z
54  After th_sph rotation, ph_sph rotation is about the original Z axis
55  In Lab coordinate system:
56
57  X = r*sin(th_geo)*cos(ph_geo) X = r*sin(th_sph)*cos(ph_sph)
58  Y = r*sin(ph_geo) Y = r*sin(th_sph)*sin(ph_sph)
59  Z = r*cos(th_geo)*cos(ph_geo) Z = r*cos(th_sph)
60
61  cos(th_sph) = cos(th_geo)*cos(ph_geo)
62  cos(ph_sph) = sin(th_geo)*cos(ph_geo)/sqrt(1-cos^2(th_geo)*cos^2(ph_geo))
63
64  GeoToSph is coded so that
65  1. negative th_geo and ph_geo = 0 returns th_sph=abs(th_geo) and ph_sph =180
66  2. positive th_geo and ph_geo = 0 returns th_sph=th_geo and ph_sph =0
67
68  Using the spherical angles, the TRotation fToLabRot and inverse fToTraRot are calculated
69  fToLabRot is rotation matrix from the spectrometer TRANSPORT system to Lab system
70  TRANSPORT coordinates are +X_tra points vertically down, +Z_tra is along the central ray and +Y_tra = ZxX
71
72  For ph_sph = 0 X_lab = Y_tra*cos(th_sph) + Z_tra*sin(th_sph)
73  Y_lab = -X_tra
74  Z_lab = -Y_tra*sin(th_sph) + Z_tra*cos(th_sph)
75  For ph_sph = 180 X_lab = Y_tra*cos(th_sph) - Z_tra*sin(th_sph)
76  Y_lab = -X_tra
77  Z_lab = Y_tra*sin(th_sph) + Z_tra*cos(th_sph)
78
79
80 */
82
83 #include "THcHallCSpectrometer.h"
84 #include "THaTrackingDetector.h"
85 #include "THcGlobals.h"
86 #include "THcParmList.h"
87 #include "THaTrack.h"
88 #include "THaTrackProj.h"
89 #include "THaTriggerTime.h"
90 #include "TMath.h"
91 #include "TList.h"
92
93 #include "THcRawShowerHit.h"
94 #include "THcSignalHit.h"
95 #include "THcShower.h"
96 #include "THcHitList.h"
97 #include "THcHodoscope.h"
98
99 #include <vector>
100 #include <cstring>
101 #include <cstdio>
102 #include <cstdlib>
103 #include <iostream>
104 #include <fstream>
105
106 using std::vector;
107 using namespace std;
108
109 //_____________________________________________________________________________
110 THcHallCSpectrometer::THcHallCSpectrometer( const char* name, const char* description ) :
111  THaSpectrometer( name, description ), fPresent(kTRUE)
112 {
113  // Constructor. Defines the standard detectors for the HRS.
114  // AddDetector( new THaTriggerTime("trg","Trigger-based time offset"));
115
116  //sc_ref = static_cast<THaScintillator*>(GetDetector("s1"));
117
119  eventtypes.clear();
120 }
121
122 //_____________________________________________________________________________
124 {
125  // Destructor
126
127  DefineVariables( kDelete );
128 }
129
130 //_____________________________________________________________________________
132 {
133  // Define/delete standard variables for a spectrometer (tracks etc.)
134  // Can be overridden or extended by derived (actual) apparatuses
135  if( mode == kDefine && fIsSetup ) return kOK;
136  THaSpectrometer::DefineVariables( mode );
137  fIsSetup = ( mode == kDefine );
138  RVarDef vars[] = {
139  { "tr.betachisq", "Chi2 of beta", "fTracks.THaTrack.GetBetaChi2()"},
140  { "tr.GoodPlane4", "Flag for track hitting hodo plane 4", "fTracks.THaTrack.GetGoodPlane4()"},
141  { "tr.GoodPlane3", "Flag for track hitting hodo plane 3", "fTracks.THaTrack.GetGoodPlane3()"},
142  { "tr.fptime", "Track hodo focal plane time", "fTracks.THaTrack.GetFPTime()"},
143  { "tr.npmt", "Track number of hodo PMTs hit", "fTracks.THaTrack.GetNPMT()"},
144  { "tr.PruneSelect", "Prune Select ID", "fPruneSelect"},
145  { "present", "Trigger Type includes this spectrometer", "fPresent"},
146  { 0 }
147  };
148
149
150  return DefineVarsFromList( vars, mode );
151 }
152
153 //_____________________________________________________________________________
155 {
156  if( set )
157  fProperties |= kSortTracks;
158  else
159  fProperties &= ~kSortTracks;
160
161  return set;
162 }
163
164 //_____________________________________________________________________________
166 {
167  return ((fProperties & kSortTracks) != 0);
168 }
169
170 //_____________________________________________________________________________
172 {
173  fNReconTerms = 0;
174  fReconTerms.clear();
175  fAngSlope_x = 0.0;
176  fAngSlope_y = 0.0;
177  fAngOffset_x = 0.0;
178  fAngOffset_y = 0.0;
179  fDetOffset_x = 0.0;
180  fDetOffset_y = 0.0;
181  fZTrueFocus = 0.0;
182 }
183 //_____________________________________________________________________________
185 {
200  static const char* const here = "THcHallCSpectrometer::ReadDatabase";
201
202 #ifdef WITH_DEBUG
203  cout << "In THcHallCSpectrometer::ReadDatabase()" << endl;
204 #endif
205
206  const char* detector_name = "hod";
207  //THaApparatus* app = GetDetector();
209  // THaDetector* det = app->GetDetector( shower_detector_name );
210
211  if( dynamic_cast<THcHodoscope*>(det) ) {
212  fHodo = static_cast<THcHodoscope*>(det); // fHodo is a membervariable
213  } else {
214  Error("THcHallCSpectrometer", "Cannot find hodoscope detector %s",
215  detector_name );
216  fHodo = NULL;
217  }
218
219
221  if( dynamic_cast<THcDC*>(detdc) ) {
222  fDC = static_cast<THcDC*>(detdc); // fHodo is a membervariable
223  } else {
224  Error("THcHallCSpectrometer", "Cannot find detector DC");
225  fDC = NULL;
226  }
227
228
229  // Get the matrix element filename from the variable store
230  // Read in the matrix
231
233
234  char prefix[2];
235
236 #ifdef WITH_DEBUG
237  cout << " GetName() " << GetName() << endl;
238 #endif
239
240  prefix[0]=tolower(GetName()[0]);
241  prefix[1]='\0';
242
243  string reconCoeffFilename;
244  DBRequest list[]={
245  {"_recon_coeff_filename", &reconCoeffFilename, kString },
246  {"theta_offset", &fThetaOffset, kDouble },
247  {"phi_offset", &fPhiOffset, kDouble },
248  {"delta_offset", &fDeltaOffset, kDouble },
249  {"thetacentral_offset", &fThetaCentralOffset, kDouble },
250  {"_oopcentral_offset", &fOopCentralOffset, kDouble },
251  {"pcentral_offset", &fPCentralOffset, kDouble },
252  {"pcentral", &fPcentral, kDouble },
253  {"satcorr", &fSatCorr, kDouble, 0, 1},
254  {"theta_lab", &fTheta_lab, kDouble },
255  {"partmass", &fPartMass, kDouble },
256  {"phi_lab", &fPhi_lab, kDouble, 0, 1},
257  {"mispointing_x", &fMispointing_x, kDouble, 0, 1},
258  {"mispointing_y", &fMispointing_y, kDouble, 0, 1},
259  {"sel_using_scin", &fSelUsingScin, kInt, 0, 1},
260  {"sel_using_prune", &fSelUsingPrune, kInt, 0, 1},
261  {"sel_ndegreesmin", &fSelNDegreesMin, kDouble, 0, 1},
262  {"sel_dedx1min", &fSeldEdX1Min, kDouble, 0, 1},
263  {"sel_dedx1max", &fSeldEdX1Max, kDouble, 0, 1},
264  {"sel_betamin", &fSelBetaMin, kDouble, 0, 1},
265  {"sel_betamax", &fSelBetaMax, kDouble, 0, 1},
266  {"sel_etmin", &fSelEtMin, kDouble, 0, 1},
267  {"sel_etmax", &fSelEtMax, kDouble, 0, 1},
268  {"hodo_num_planes", &fNPlanes, kInt },
269  {"scin_2x_zpos", &fScin2XZpos, kDouble, 0, 1},
270  {"scin_2x_dzpos", &fScin2XdZpos, kDouble, 0, 1},
271  {"scin_2y_zpos", &fScin2YZpos, kDouble, 0, 1},
272  {"scin_2y_dzpos", &fScin2YdZpos, kDouble, 0, 1},
273  {"prune_xp", &fPruneXp, kDouble, 0, 1},
274  {"prune_yp", &fPruneYp, kDouble, 0, 1},
275  {"prune_ytar", &fPruneYtar, kDouble, 0, 1},
276  {"prune_delta", &fPruneDelta, kDouble, 0, 1},
277  {"prune_beta", &fPruneBeta, kDouble, 0, 1},
278  {"prune_df", &fPruneDf, kDouble, 0, 1},
279  {"prune_chibeta", &fPruneChiBeta, kDouble, 0, 1},
280  {"prune_npmt", &fPruneNPMT, kDouble, 0, 1},
281  {"prune_fptime", &fPruneFpTime, kDouble, 0, 1},
282  {"prune_DipoleExit", &fPruneDipoleExit, kDouble, 0, 1},
283  {0}
284  };
285
286  // Default values
288  fSelUsingScin = 0;
289  fSelUsingPrune = 0;
290  fPruneXp = .2;
291  fPruneYp = .2;
292  fPruneYtar = 20.;
293  fPruneDelta = 30.;
294  fPruneBeta = 30.;
295  fPruneDf= 1;
296  fPruneChiBeta= 100.;
297  fPruneNPMT= 6;
298  fPruneFpTime= 1000.;
299  fPhi_lab = 0.;
300  fSatCorr=0.;
301  fMispointing_x=999.;
302  fMispointing_y=999.;
306  if (prefix[0]=='h') fUseHMSDipoleExitWindow=kTRUE;
307  if (prefix[0]=='p') fUseSHMSDipoleExitWindow=kTRUE;
308
309  // mispointing in transport system y is horizontal and +x is vertical down
310  if (fMispointing_y == 999.) {
311
312  if (prefix[0]=='h') {
313
315  else {fMispointing_y = 0.1*(0.52-0.012*40.+0.002*40.*40.);}
316
317  gHcParms->Define("hmispointing_y","HMS Y-Mispointing", fMispointing_y);
318  }
319
320  if (prefix[0]=='p') {
321
322  fMispointing_y = 0.1*(-0.6);
323  gHcParms->Define("pmispointing_y","SHMS Y-Mispointing", fMispointing_y);
324  cout << prefix[0] << " From Formula Mispointing_y = " << fMispointing_y << endl;
325  }
326
327  }
328
329  else {
330  cout << prefix[0] << " From Parameter Set Mispointing_y = " << fMispointing_y << endl;
331  }
332
333
334  if (fMispointing_x == 999.) {
335
336  if (prefix[0]=='h') {
337
338
340  else {fMispointing_x = 0.1*(2.37-0.086*50+0.0012*50.*50.);}
341
342  gHcParms->Define("hmispointing_x","HMS X-Mispointing", fMispointing_x);
343  cout << prefix[0] << " From Formula Mispointing_x = " << fMispointing_x << endl;
344
345  }
346
347  if (prefix[0]=='p') {
348
349  fMispointing_x = 0.1*(-1.26);
350
351  gHcParms->Define("pmispointing_x","SHMS X-Mispointing", fMispointing_x);
352  cout << prefix[0] << " From Formula Mispointing_x = " << fMispointing_x << endl;
353
354  }
355
356  }
357
358  else {
359  cout << prefix[0] << " From Parameter Set Mispointing_x = " << fMispointing_x << endl;
360  }
361  //
362
363
364  //EnforcePruneLimits();
365
366 #ifdef WITH_DEBUG
367  cout << "\n\n\nhodo planes = " << fNPlanes << endl;
368  cout << "sel using scin = " << fSelUsingScin << endl;
369  cout << "fPruneXp = " << fPruneXp << endl;
370  cout << "fPruneYp = " << fPruneYp << endl;
371  cout << "fPruneYtar = " << fPruneYtar << endl;
372  cout << "fPruneDelta = " << fPruneDelta << endl;
373  cout << "fPruneBeta = " << fPruneBeta << endl;
374  cout << "fPruneDf = " << fPruneDf << endl;
375  cout << "fPruneChiBeta = " << fPruneChiBeta << endl;
376  cout << "fPruneFpTime = " << fPruneFpTime << endl;
377  cout << "fPruneNPMT = " << fPruneNPMT << endl;
378  cout << "sel using prune = " << fSelUsingPrune << endl;
379 #endif
380  cout << "fPartMass = " << fPartMass << endl;
381  cout << "fPcentral = " << fPcentral << " " <<fPCentralOffset << endl;
382  cout << "fThetalab = " << fTheta_lab << " " <<fThetaCentralOffset << endl;
383  fPcentral= fPcentral*(1.+fPCentralOffset/100.);
384  // Check that these offsets are in radians
387  // SetCentralAngles method in podd THaSpectrometer
388  // fTheta_lab and ph are geographical angles, converts to spherical coordinates
389  // Need to set fTheta_lab to negative for spectrometer like HMS on beam right
390  // This gives phi_sph = 180 th_sph=abs(th_geo)
391  // Computes TRotation fToLabRot and fToTraRot
392  Bool_t bend_down = kFALSE;
393  SetCentralAngles(fTheta_lab, ph, bend_down);
394  Double_t off_z = 0.0;
395  fPointingOffset.SetXYZ( fMispointing_x, fMispointing_y, off_z );
396  //
397  ifstream ifile;
398  ifile.open(reconCoeffFilename.c_str());
399  if(!ifile.is_open()) {
400  Error(here, "error opening reconstruction coefficient file %s",reconCoeffFilename.c_str());
401  // return kInitError; // Is this the right return code?
402  return kOK;
403  }
404
405  string line="!";
406  int good=1;
407  while(good && line[0]=='!') {
408  good = getline(ifile,line).good();
409  // cout << line << endl;
410  }
411  // Read in focal plane rotation coefficients
412  // Probably not used, so for now, just paste in fortran code as a comment
413  while(good && line.compare(0,4," ---")!=0) {
421  good = getline(ifile,line).good();
422  }
423  // Read in reconstruction coefficients and exponents
424  line=" ";
425  good = getline(ifile,line).good();
426  // cout << line << endl;
427  fNReconTerms = 0;
428  fReconTerms.clear();
429  fReconTerms.reserve(500);
430  //cout << "Reading matrix elements" << endl;
431  while(good && line.compare(0,4," ---")!=0) {
432  fReconTerms.push_back(reconTerm());
433  sscanf(line.c_str()," %le %le %le %le %1d%1d%1d%1d%1d"
434  ,&fReconTerms[fNReconTerms].Coeff[0],&fReconTerms[fNReconTerms].Coeff[1]
435  ,&fReconTerms[fNReconTerms].Coeff[2],&fReconTerms[fNReconTerms].Coeff[3]
436  ,&fReconTerms[fNReconTerms].Exp[0]
437  ,&fReconTerms[fNReconTerms].Exp[1]
438  ,&fReconTerms[fNReconTerms].Exp[2]
439  ,&fReconTerms[fNReconTerms].Exp[3]
440  ,&fReconTerms[fNReconTerms].Exp[4]);
441  fNReconTerms++;
442  good = getline(ifile,line).good();
443  }
444  cout << "Read " << fNReconTerms << " matrix element terms" << endl;
445  if(!good) {
446  Error(here, "Error processing reconstruction coefficient file %s",reconCoeffFilename.c_str());
447  return kInitError; // Is this the right return code?
448  }
449  return kOK;
450 }
451
452 //_____________________________________________________________________________
454 {
455
456  fPruneXp = TMath::Max( 0.08, fPruneXp);
457  fPruneYp = TMath::Max( 0.04, fPruneYp);
461  fPruneDf = TMath::Max( 1.0, fPruneDf);
465 }
466
467 //_____________________________________________________________________________
469 {
482  fNtracks = tracks.GetLast()+1;
483
484  /*
485  for (Int_t it=0;it<tracks.GetLast()+1;it++) {
486  THaTrack* track = static_cast<THaTrack*>( tracks[it] );
487  Double_t xptar=kBig,yptar=kBig,ytar=kBig,delta=kBig;
488  Double_t xtar=0;
489  CalculateTargetQuantities(track,xtar,xptar,ytar,yptar,delta);
490  // Transfer results to track
491  // No beam raster yet
492  //; In transport coordinates phi = hyptar = dy/dz and theta = hxptar = dx/dz
493  //; but for unknown reasons the yp offset is named htheta_offset
494  //; and the xp offset is named hphi_offset
495
496  track->SetTarget(0.0, ytar*100.0, xptar, yptar);
497  track->SetDp(delta*100.0); // Percent.
498  Double_t ptemp = fPcentral*(1+track->GetDp()/100.0);
499  cout << " it = " << " Set track p = " << ptemp << endl;
500  track->SetMomentum(ptemp);
501  TVector3 pvect_temp;
502  TransportToLab(track->GetP(),track->GetTTheta(),track->GetTPhi(),pvect_temp);
503  track->SetPvect(pvect_temp);
504  }
505  */
506  fPruneSelect=-1.;
507  if (fHodo==0 || (( fSelUsingScin == 0 ) && ( fSelUsingPrune == 0 )) ) {
508  BestTrackSimple();
509  } else if (fHodo!=0 && fSelUsingPrune !=0) {
511  } else if (fHodo!=0){
513  }
514
515
516  return 0;
517 }
518 //
519 void THcHallCSpectrometer::CalculateTargetQuantities(THaTrack* track,Double_t& xtar,Double_t& xptar,Double_t& ytar,Double_t& yptar,Double_t& delta)
520 {
529  Double_t hut[5];
530  Double_t hut_rot[5];
531
532  hut[0] = track->GetX()/100.0 + fZTrueFocus*track->GetTheta() + fDetOffset_x;//m
533  hut[1] = track->GetTheta() + fAngOffset_x;//radians
534  hut[2] = track->GetY()/100.0 + fZTrueFocus*track->GetPhi() + fDetOffset_y;//m
535  hut[3] = track->GetPhi() + fAngOffset_y;//radians
536
537  hut[4] = xtar/100.0;
538
539  // Retrieve the focal plane coordnates
540  // Do the transformation
541  // Stuff results into track
542  hut_rot[0] = hut[0];
543  hut_rot[1] = hut[1] + hut[0]*fAngSlope_x;
544  hut_rot[2] = hut[2];
545  hut_rot[3] = hut[3] + hut[2]*fAngSlope_y;
546  hut_rot[4] = hut[4];
547
548  // Compute COSY sums
549  Double_t sum[4];
550  for(Int_t k=0;k<4;k++) {
551  sum[k] = 0.0;
552  }
553  for(Int_t iterm=0;iterm<fNReconTerms;iterm++) {
554  Double_t term=1.0;
555  for(Int_t j=0;j<5;j++) {
556  if(fReconTerms[iterm].Exp[j]!=0) {
557  term *= pow(hut_rot[j],fReconTerms[iterm].Exp[j]);
558  }
559  }
560  for(Int_t k=0;k<4;k++) {
561  sum[k] += term*fReconTerms[iterm].Coeff[k];
562  }
563  }
564  xptar=sum[0] + fPhiOffset;
565  ytar=sum[1];
566  yptar=sum[2] + fThetaOffset;
567  delta=sum[3] + fDeltaOffset;
568  if (fSatCorr == 2000) {
569  Double_t p0corr = 0.82825*fPcentral-1.223 ;
570  delta = delta + p0corr*xptar/100.;
571  }
572 }
573 //
574 //_____________________________________________________________________________
576 {
577  if( fNtracks > 0 ) {
578  Int_t hit_gold_track=0; // find track with index =0 which is best track
579  Int_t hit_dc_track=1; //
580  for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){
581  THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) );
582  if (aTrack->GetIndex()==0) {
583  hit_gold_track=itrack;
584  hit_dc_track = aTrack->GetTrkNum();
585  }
586  }
587  fDC->SetFocalPlaneBestTrack(hit_dc_track-1);
588  fGoldenTrack = static_cast<THaTrack*>( fTracks->At(hit_gold_track) );
589  fTrkIfo = *fGoldenTrack;
590  fTrk = fGoldenTrack;
591  } else {
592  fGoldenTrack = NULL;
593  }
594
595  return TrackTimes( fTracks );
596 }
597
598 //_____________________________________________________________________________
600 {
606  if( GetTrSorting() ) fTracks->Sort();
607
609  for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){
610  THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) );
611  aTrack->SetIndex(1);
612  if (itrack==0) aTrack->SetIndex(0);
613  }
614
615  return(0);
616 }
617
618 //_____________________________________________________________________________
620 {
628  Double_t chi2Min;
629
630  if( fNtracks > 0 ) {
631  fGoodTrack = -1;
632  chi2Min = 10000000000.0;
633  Int_t y2Dmin = 100;
634  Int_t x2Dmin = 100;
635
636  Bool_t* x2Hits = new Bool_t [fHodo->GetNPaddles(2)];
637  Bool_t* y2Hits = new Bool_t [fHodo->GetNPaddles(3)];
638  for (UInt_t j = 0; j < fHodo->GetNPaddles(2); j++ ){
639  x2Hits[j] = kFALSE;
640  }
641  for (UInt_t j = 0; j < fHodo->GetNPaddles(3); j++ ){
642  y2Hits[j] = kFALSE;
643  }
644
645  for (Int_t rawindex=0; rawindex<fHodo->GetTotHits(); rawindex++) {
646  Int_t ip = fHodo->GetGoodRawPlane(rawindex);
647  if(ip==2) { // X2
650  } else if (ip==3) { // Y2
653  }
654  }
655
656  for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){
657  Double_t chi2PerDeg;
658
659  THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) );
660  if (!aTrack) {delete[] x2Hits; delete[] y2Hits; return -1;}
661
662  if ( aTrack->GetNDoF() > fSelNDegreesMin ){
663  chi2PerDeg = aTrack->GetChi2() / aTrack->GetNDoF();
664
665  if( ( aTrack->GetDedx() > fSeldEdX1Min ) &&
666  ( aTrack->GetDedx() < fSeldEdX1Max ) &&
667  ( aTrack->GetBeta() > fSelBetaMin ) &&
668  ( aTrack->GetBeta() < fSelBetaMax ) &&
669  ( aTrack->GetEnergy() > fSelEtMin ) &&
670  ( aTrack->GetEnergy() < fSelEtMax ) ) {
671  Int_t x2D, y2D;
672
673  if ( fNtracks > 1 ){
674  Double_t hitpos3 = aTrack->GetX() + aTrack->GetTheta() * ( fScin2XZpos + 0.5 * fScin2XdZpos );
675  Int_t icounter3 = TMath::Nint( ( hitpos3 - fHodo->GetPlaneCenter(2) ) / fHodo->GetPlaneSpacing(2) ) + 1;
676  Int_t hitCnt3 = TMath::Max( TMath::Min(icounter3, (Int_t) fHodo->GetNPaddles(2) ) , 1); // scin_2x_nr = 16
677  // fHitDist3 = fHitPos3 - ( fHodo->GetPlaneSpacing(2) * ( hitCnt3 - 1 ) + fHodo->GetPlaneCenter(2) );
678  Double_t hitpos4 = aTrack->GetY() + aTrack->GetPhi() * ( fScin2YZpos + 0.5 * fScin2YdZpos );
679  Int_t icounter4 = TMath::Nint( ( fHodo->GetPlaneCenter(3) - hitpos4 ) / fHodo->GetPlaneSpacing(3) ) + 1;
680  Int_t hitCnt4 = TMath::Max( TMath::Min(icounter4, (Int_t) fHodo->GetNPaddles(3) ) , 1); // scin_2y_nr = 10
681  // fHitDist4 = fHitPos4 - ( fHodo->GetPlaneCenter(3) - fHodo->GetPlaneSpacing(3) * ( hitCnt4 - 1 ) );
682  // Plane 3
683  Int_t mindiff=1000;
684  for (UInt_t i = 0; i < fHodo->GetNPaddles(2); i++ ){
685  if ( x2Hits[i] ) {
686  Int_t diff = TMath::Abs((Int_t)hitCnt3-(Int_t)i-1);
687  if (diff < mindiff) mindiff = diff;
688  }
689  }
690  if(mindiff < 1000) {
691  x2D = mindiff;
692  } else {
693  x2D = 0; // Is this what we really want if there were no hits on this plane?
694  }
695
696  // Plane 4
697  mindiff = 1000;
698  for (UInt_t i = 0; i < fHodo->GetNPaddles(3); i++ ){
699  if ( y2Hits[i] ) {
700  Int_t diff = TMath::Abs((Int_t)hitCnt4-(Int_t)i-1);
701  if (diff < mindiff) mindiff = diff;
702  }
703  }
704  if(mindiff < 1000) {
705  y2D = mindiff;
706  } else {
707  y2D = 0; // Is this what we really want if there were no hits on this plane?
708  }
709  } else { // Only a single track
710  x2D = 0.;
711  y2D = 0.;
712  }
713
714  if ( y2D <= y2Dmin ) {
715  if ( y2D < y2Dmin ) {
716  x2Dmin = 100;
717  chi2Min = 10000000000.;
718  } // y2D min
719
720  if ( x2D <= x2Dmin ){
721  if ( x2D < x2Dmin ){
722  chi2Min = 10000000000.0;
723  } // condition x2D
724  if ( chi2PerDeg < chi2Min ){
725
726  fGoodTrack = itrack; // fGoodTrack = itrack
727  y2Dmin = y2D;
728  x2Dmin = x2D;
729  chi2Min = chi2PerDeg;
730
731  fGoldenTrack = static_cast<THaTrack*>( fTracks->At( fGoodTrack ) );
732  fTrkIfo = *fGoldenTrack;
733  fTrk = fGoldenTrack;
734  }
735  } // condition x2D
736  } // condition y2D
737  } // conditions for dedx, beta and trac energy
738  } // confition for fNFreeFP greater than fSelNDegreesMin
739  } // loop over tracks
740
741  // Fall back to using track with minimum chi2
742  if ( fGoodTrack == -1 ){
743
744  chi2Min = 10000000000.0;
745  for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){
746  Double_t chi2PerDeg;
747  THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) );
748  if (!aTrack) return -1;
749
750  if ( aTrack->GetNDoF() > fSelNDegreesMin ){
751  chi2PerDeg = aTrack->GetChi2() / aTrack->GetNDoF();
752
753  if ( chi2PerDeg < chi2Min ){
754
755  fGoodTrack = iitrack;
756  chi2Min = chi2PerDeg;
757
758  }
759  }
760  } // loop over trakcs
761  // Set index for fGoodTrack = 0
762  for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){
763  THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) );
764  aTrack->SetIndex(1);
765  if (iitrack==fGoodTrack) aTrack->SetIndex(0);
766  }
767  //
768  }
769
770  }
771
772  return(0);
773 }
774
775 //_____________________________________________________________________________
777 {
786  Int_t nGood;
787  Double_t chi2Min;
788
789  if ( fNtracks > 0 ) {
790  chi2Min = 10000000000.0;
791  fGoodTrack = 0;
792  vector<bool> keep(fNtracks);
793  vector<Int_t> reject(fNtracks);
794
795  THaTrack *testTracks[fNtracks];
796
797  // ! Initialize all tracks to be good
798  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
799  keep[ptrack] = kTRUE;
800  reject[ptrack] = 0;
801  testTracks[ptrack] = static_cast<THaTrack*>( fTracks->At(ptrack) );
802  if (!testTracks[ptrack]) return -1;
803  }
804  fPruneSelect = 0;
805  Double_t PruneSelect=0;
806  // ! Prune on xptar
807  nGood = 0;
808  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
809  if ( ( TMath::Abs( testTracks[ptrack]->GetTTheta() ) < fPruneXp ) && ( keep[ptrack] ) ){
810  nGood ++;
811  }
812  }
813  if ( nGood > 0 ) {
814  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
815  if ( TMath::Abs( testTracks[ptrack]->GetTTheta() ) >= fPruneXp ){
816  keep[ptrack] = kFALSE;
817  reject[ptrack] += 1;
818  }
819  }
820  }
821  PruneSelect++;
822  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
823  // ! Prune on yptar
824  nGood = 0;
825  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
826  if ( ( TMath::Abs( testTracks[ptrack]->GetTPhi() ) < fPruneYp ) && ( keep[ptrack] ) ){
827  nGood ++;
828  }
829  }
830  if (nGood > 0 ) {
831  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
832  if ( TMath::Abs( testTracks[ptrack]->GetTPhi() ) >= fPruneYp ){
833  keep[ptrack] = kFALSE;
834  reject[ptrack] += 2;
835
836  }
837  }
838  }
839  PruneSelect++;
840  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
841
842  // ! Prune on ytar
843  nGood = 0;
844  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
845  if ( ( TMath::Abs( testTracks[ptrack]->GetTY() ) < fPruneYtar ) && ( keep[ptrack] ) ){
846  nGood ++;
847  }
848  }
849  if (nGood > 0 ) {
850  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
851  if ( TMath::Abs( testTracks[ptrack]->GetTY() ) >= fPruneYtar ){
852  keep[ptrack] = kFALSE;
853  reject[ptrack] += 10;
854  }
855  }
856  }
857  PruneSelect++;
858  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
859
860  // ! Prune on delta
861  nGood = 0;
862  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
863  if ( ( TMath::Abs( testTracks[ptrack]->GetDp() ) < fPruneDelta ) && ( keep[ptrack] ) ){
864  nGood ++;
865  }
866  }
867  if (nGood > 0 ) {
868  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
869  if ( TMath::Abs( testTracks[ptrack]->GetDp() ) >= fPruneDelta ){
870  keep[ptrack] = kFALSE;
871  reject[ptrack] += 20;
872  }
873  }
874  }
875  PruneSelect++;
876  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
877
878  // ! Prune on dipole exit
879  nGood = 0;
880  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
881  Double_t xfp=testTracks[ptrack]->GetX();
882  Double_t yfp=testTracks[ptrack]->GetY();
883  Double_t xpfp=testTracks[ptrack]->GetTheta();
884  Double_t ypfp=testTracks[ptrack]->GetPhi();
885  if ( fPruneDipoleExit==1 && InsideDipoleExitWindow(xfp,xpfp,yfp,ypfp) && ( keep[ptrack] ) ){
886  nGood ++;
887  }
888  }
889  if (nGood > 0 ) {
890  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
891  Double_t xfp=testTracks[ptrack]->GetX();
892  Double_t yfp=testTracks[ptrack]->GetY();
893  Double_t xpfp=testTracks[ptrack]->GetTheta();
894  Double_t ypfp=testTracks[ptrack]->GetPhi();
895  if (!InsideDipoleExitWindow(xfp,xpfp,yfp,ypfp) ){
896  keep[ptrack] = kFALSE;
897  reject[ptrack] += 30;
898  }
899  }
900  }
901  PruneSelect++;
902  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
903
904  // ! Prune on beta
905  nGood = 0;
906  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
907  Double_t p = testTracks[ptrack]->GetP();
908  Double_t betaP = p / TMath::Sqrt( p * p + fPartMass * fPartMass );
909  if ( ( TMath::Abs( testTracks[ptrack]->GetBeta() - betaP ) < fPruneBeta ) && ( keep[ptrack] ) ){
910  nGood ++;
911  }
912  }
913  if (nGood > 0 ) {
914  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
915  Double_t p = testTracks[ptrack]->GetP();
916  Double_t betaP = p / TMath::Sqrt( p * p + fPartMass * fPartMass );
917  if ( TMath::Abs( testTracks[ptrack]->GetBeta() - betaP ) >= fPruneBeta ) {
918  keep[ptrack] = kFALSE;
919  reject[ptrack] += 100;
920  }
921  }
922  }
923  PruneSelect++;
924  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
925
926  // ! Prune on deg. freedom for track chisq
927  nGood = 0;
928  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
929  if ( ( testTracks[ptrack]->GetNDoF() >= fPruneDf ) && ( keep[ptrack] ) ){
930  nGood ++;
931  }
932  }
933  if (nGood > 0 ) {
934  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
935  if ( testTracks[ptrack]->GetNDoF() < fPruneDf ){
936  keep[ptrack] = kFALSE;
937  reject[ptrack] += 200;
938  }
939  }
940  }
941
942  PruneSelect++;
943  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
945  nGood = 0;
946  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
947  if ( ( testTracks[ptrack]->GetNPMT() >= fPruneNPMT ) && ( keep[ptrack] ) ){
948  nGood ++;
949  }
950  }
951  if (nGood > 0 ) {
952  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
953  if ( testTracks[ptrack]->GetNPMT() < fPruneNPMT ){
954  keep[ptrack] = kFALSE;
955  reject[ptrack] += 100000;
956  }
957  }
958  }
959  PruneSelect++;
960  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
961
962  // ! Prune on beta chisqr
963  nGood = 0;
964  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
965  if ( ( testTracks[ptrack]->GetBetaChi2() < fPruneChiBeta ) &&
966  ( testTracks[ptrack]->GetBetaChi2() > 0.01 ) && ( keep[ptrack] ) ){
967  nGood ++;
968  }
969  }
970  if (nGood > 0 ) {
971  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
972  if ( ( testTracks[ptrack]->GetBetaChi2() >= fPruneChiBeta ) ||
973  ( testTracks[ptrack]->GetBetaChi2() <= 0.01 ) ){
974  keep[ptrack] = kFALSE;
975  reject[ptrack] += 1000;
976  }
977  }
978  }
979  PruneSelect++;
980  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
981
982  // ! Prune on fptime
983  nGood = 0;
984  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
985  if ( ( TMath::Abs( testTracks[ptrack]->GetFPTime() - fHodo->GetStartTimeCenter() ) < fPruneFpTime ) &&
986  ( keep[ptrack] ) ){
987  nGood ++;
988  }
989  }
990  if (nGood > 0 ) {
991  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
992  if ( TMath::Abs( testTracks[ptrack]->GetFPTime() - fHodo->GetStartTimeCenter() ) >= fPruneFpTime ) {
993  keep[ptrack] = kFALSE;
994  reject[ptrack] += 2000;
995  }
996  }
997  }
998  PruneSelect++;
999  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
1000
1001  // ! Prune on Y2 being hit
1002  nGood = 0;
1003  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
1004  if ( ( testTracks[ptrack]->GetGoodPlane4() == 1 ) && keep[ptrack] ) {
1005  nGood ++;
1006  }
1007  }
1008  if (nGood > 0 ) {
1009  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
1010  if ( testTracks[ptrack]->GetGoodPlane4() != 1 ) {
1011  keep[ptrack] = kFALSE;
1012  reject[ptrack] += 10000;
1013  }
1014  }
1015  }
1016  PruneSelect++;
1017  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
1018
1019  // ! Prune on X2 being hit
1020  nGood = 0;
1021  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
1022  if ( ( testTracks[ptrack]->GetGoodPlane3() == 1 ) && keep[ptrack] ) {
1023  nGood ++;
1024  }
1025  }
1026  if (nGood > 0 ) {
1027  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
1028  if ( testTracks[ptrack]->GetGoodPlane3() != 1 ) {
1029  keep[ptrack] = kFALSE;
1030  reject[ptrack] += 20000;
1031  }
1032  }
1033  }
1034  PruneSelect++;
1035  if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
1036
1037  // ! Pick track with best chisq if more than one track passed prune tests
1038  Double_t chi2PerDeg = 0.;
1039  for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
1040
1041  chi2PerDeg = testTracks[ptrack]->GetChi2() / testTracks[ptrack]->GetNDoF();
1042
1043  if ( ( chi2PerDeg < chi2Min ) && ( keep[ptrack] ) ){
1044  fGoodTrack = ptrack;
1045  chi2Min = chi2PerDeg;
1046  }
1047  }
1048  PruneSelect++;
1049  if (fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
1050  // Set index=0 for fGoodTrack
1051  for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){
1052  THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) );
1053  aTrack->SetIndex(1);
1054  if (iitrack==fGoodTrack) aTrack->SetIndex(0);
1055  }
1056  //
1057
1058
1059  }
1060
1061  return(0);
1062 }
1063
1064 //_____________________________________________________________________________
1066  // Do the actual track-timing (beta) calculation.
1067  // Use multiple scintillators to average together and get "best" time at S1.
1068  //
1069  // To be useful, a meaningful timing resolution should be assigned
1070  // to each Scintillator object (part of the database).
1071
1072
1073  return 0;
1074 }
1075
1077 {
1078
1079  fPresent=kTRUE;
1080  if(eventtypes.size()!=0) {
1081  Int_t evtype = evdata.GetEvType();
1082  if(!IsMyEvent(evtype)) {
1083  fPresent = kFALSE;
1084  }
1085  }
1086
1087  return THaSpectrometer::Decode(evdata);
1088 }
1089
1090 //_____________________________________________________________________________
1092 {
1093  // Override THaSpectrometer with nop method. All needed kinamatics
1095
1096  return kOK;
1097 }
1098
1100  eventtypes.push_back(evtype);
1101 }
1102
1104  eventtypes.clear();
1106 }
1107
1109 {
1110  for (UInt_t i=0; i < eventtypes.size(); i++) {
1111  if (evtype == eventtypes[i]) return kTRUE;
1112  }
1113
1114  return kFALSE;
1115 }
1116 //
1118  Bool_t inside=kTRUE;
1119  Double_t DipoleExitWindowZpos=0.;
1120  if (fUseSHMSDipoleExitWindow) DipoleExitWindowZpos=-307.;
1121  if (fUseHMSDipoleExitWindow) DipoleExitWindowZpos=-147.48;
1122  Double_t xdip = x_fp + xp_fp*DipoleExitWindowZpos;
1123  Double_t ydip = y_fp + yp_fp*DipoleExitWindowZpos;
1124  if (fUseSHMSDipoleExitWindow) inside = SHMSDipoleExitWindow(xdip,ydip);
1125  if (fUseHMSDipoleExitWindow) inside = HMSDipoleExitWindow(xdip,ydip);
1126  return inside;
1127 }
1128 //
1130  Bool_t insideSHMS=kTRUE;
1133  Double_t hwid=11.549/2.;
1134  if ( TMath::Abs(ydip) < hwid) {
1135  if (TMath::Abs(xdip) > (crad+voffset)) insideSHMS=kFALSE;
1136  } else {
1137  if ( ydip >=hwid) {
1139  }
1140  if ( ydip <=-hwid) {
1142  }
1143  }
1144  return insideSHMS;
1145  }
1146 //
1148  Bool_t insideHMS=kTRUE;
1149  Double_t xpipe_offset = 2.8;
1150  Double_t ypipe_offset = 0.0;
1153  return insideHMS ;
1154 }
1155
1156 //_____________________________________________________________________________
Loads parameters to characterize a Hall C spectrometer.
std::vector< reconTerm > fReconTerms
virtual Int_t BestTrackUsingScin()
Choose best track using closeness to scintillator hits.
TLine * line
Int_t GetLast() const
Bool_t SHMSDipoleExitWindow(Double_t x_dip, Double_t y_dip)
Double_t GetPlaneCenter(Int_t ip)
Definition: THcHodoscope.h:114
virtual Int_t DefineVariables(EMode mode=kDefine)
virtual void CalculateTargetQuantities(THaTrack *track, Double_t &gbeam_y, Double_t &xptar, Double_t &ytar, Double_t &yptar, Double_t &delta)
Transport focal plane track to target.
Short_t Min(Short_t a, Short_t b)
int Int_t
bool Bool_t
Bool_t SetTrSorting(Bool_t set=kFALSE)
STL namespace.
std::vector< Int_t > eventtypes
Short_t Abs(Short_t d)
Analyze a package of horizontal drift chambers.
Definition: THcDC.h:23
virtual Int_t FindVertices(TClonesArray &tracks)
Reconstruct target coordinates.
virtual Int_t Decode(const THaEvData &)
Int_t GetTotHits()
Definition: THcHodoscope.h:108
double pow(double, double)
virtual void Error(const char *method, const char *msgfmt,...) const
virtual Int_t TrackTimes(TClonesArray *tracks)
void SetFocalPlaneBestTrack(Int_t golden_track_index)
Definition: THcDC.cxx:702
unsigned int UInt_t
Bool_t InsideDipoleExitWindow(Double_t x_fp, Double_t xp_fp, Double_t y_fp, Double_t yp_fp)
virtual Int_t BestTrackUsingPrune()
Choose best track after pruning.
Double_t GetPlaneSpacing(Int_t ip)
Definition: THcHodoscope.h:115
THcHallCSpectrometer(const char *name, const char *description)
static const UInt_t kSortTracks
tuple list
Definition: SConscript.py:9
const Bool_t kFALSE
virtual Int_t BestTrackSimple()
Choose best track based on Chisq.
Definition: THcHodoscope.h:105
Double_t Exp(Double_t x)
Int_t GetGoodRawPlane(Int_t iii)
Definition: THcHodoscope.h:106
virtual const char * GetName() const
Int_t LoadParmValues(const DBRequest *list, const char *prefix="")
Retrieve parameter values from the parameter cache.
double Double_t
Bool_t HMSDipoleExitWindow(Double_t x_dip, Double_t y_dip)
Definition: THcHodoscope.h:112
UInt_t GetEvType() const
ClassImp(THcDCLookupTTDConv) THcDCLookupTTDConv
Double_t GetStartTimeCenter() const
Definition: THcHodoscope.h:101
virtual void EnforcePruneLimits()
Enforce minimum values for the prune cuts.