GapWidth_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: GapWidth
3 // Module Type: analyzer
4 // File: GapWidth_module.cc
5 // Author: Tristan Blackburn (Sussex) t.blackburn@sussex.ac.uk
6 //
7 // Module skeleton is AnaTree - K. Warburton (Sheffield)
8 ////////////////////////////////////////////////////////////////////////
9 
13 #include "fhiclcpp/ParameterSet.h"
18 #include "art_root_io/TFileService.h"
19 #include "art_root_io/TFileDirectory.h"
21 
22 // LArSoft includes
34 #include "lardata/ArtDataHelper/TrackUtils.h" // lar::util::TrackPitchInView()
38 //Commenting this out since it doesn't appear to be used.
39 //#include "larsim/MCCheater/BackTrackerService.h"
40 //#include "larsim/MCCheater/ParticleInventoryService.h"
46 
47 // ROOT includes
48 #include "TTree.h"
49 #include "TTimeStamp.h"
50 #include "TLorentzVector.h"
51 #include "TH2F.h"
52 #include "TFile.h"
53 #include "TH1.h"
54 #include "TH1D.h"
55 #include "TF1.h"
56 #include "TStyle.h"
57 #include "TMath.h"
58 #include "TMinuit.h"
59 #include "TFitter.h"
60 #include "TPaveStats.h"
61 #include "TROOT.h"
62 
63 //standard library includes
64 #include <map>
65 #include <iostream>
66 #include <iomanip>
67 #include <sstream>
68 #include <cmath>
69 #include <memory>
70 #include <limits> // std::numeric_limits<>
71 #include <vector>
72 #include <algorithm>
73 #include <array>
74 #include <iterator>
75 #include <cmath>
76 
77 const int kMaxTrack = 1000; //maximum number of tracks
78 const int kMaxHits = 10000; //maximum number of hits
79 //const int kMaxClust = 10000; //maximum number of clusters // unused
80 const int kMaxTrackHits = 1000; //maximum number of space points
81 const int kMaxEvent = 10000; //maximum number of events
82 
83 //Some loop parameters
84 
85 double track1x = 0;
86 double track1y = 0;
87 double track1z = 0;
88 
89 double track2x = 0;
90 double track2y = 0;
91 double track2z = 0;
92 
93 double track1grad = 0;
94 double track2grad = 0;
95 
96 double minx = 0;
97 double miny = 0;
98 double minz = 0;
99 
100 int posminx;
101 int posminy;
103 
104 double zdiff = 0;
105 double xdiff = 0;
106 double ydiff = 0;
107 
108 double z = 0;
109 double y = 0;
110 double translation = 0;
111 
112 
113 namespace GapWidth {
114  class GapWidth;
115 
116  double xextrap(double z) {
117  if (track1grad <= 0) return fabs((track1x - (track1grad)*(minz + z)) - track2x);
118  else return fabs((track1x + (track1grad)*(minz + z)) - track2x);
119  }
120 
121  void minuitFunctionx(int& nDim, double* gout, double& result, double par[], int flg){
122  result = xextrap(par[0]);
123  }
124 
125  double zextrap(double y) {
126  if (track1grad <= 0) return fabs((track1z - (track1grad)*(miny + y)) - track2z);
127  else return fabs((track1z + (track1grad)*(miny + y)) - track2z);
128  }
129 
130  void minuitFunctionz(int& nDim, double* gout, double& result, double par[], int flg){
131  result = zextrap(par[0]);
132  }
133 }
134 
135 
137 public:
138  explicit GapWidth(fhicl::ParameterSet const & p);
139  virtual ~GapWidth();
140 
141  void analyze(art::Event const & e) override;
142  void beginJob() override;
143  void endJob() override;
144 
145 
146 private:
147 
148  void ResetVars();
149 
150  TTree* fTree;
151 
152  // Run information
153  int run;
154  int subrun;
155  int event;
156  double evttime;
157  double efield[3];
158  int t0;
159 
160  int ntracks_reco; //number of reconstructed tracks
167 
168  double trkd2[kMaxTrack];
169  double trklen[kMaxTrack];
172 
178  double trkphi[kMaxTrack];
183 
184  double trkdedx2[kMaxTrack][3][1000];
185  double trkdqdx[kMaxTrack][3][1000];
186  double trkpitchHit[kMaxTrack][3][1000];
187  double trkkinE[kMaxTrack][3];
188  double trkrange[kMaxTrack][3];
189  double trkTPC[kMaxTrack][3][1000];
190  double trkplaneid[kMaxTrack][3][1000];
191  double trkresrg[kMaxTrack][3][1000];
192  double trkPosx[kMaxTrack][3][1000];
193  double trkPosy[kMaxTrack][3][1000];
194  double trkPosz[kMaxTrack][3][1000];
195 
206 
207  double trkpitch[kMaxTrack][3];
208  int nhits;
209  int nhits2;
210  int nclust;
211 
218  double hit_ph[kMaxHits];
220 
222 
231 
232  //double fElectronsToGeV; // conversion factor // unused
233  art::ServiceHandle<geo::Geometry> fGeometry; // pointer to Geometry service
234 
239 
240  //vectors of possible minimised distances from possible start and end points between two adjacent TPC tracks
241  std::vector<double> possiblex;
242  std::vector<double> possibley;
243  std::vector<double> possiblez;
244 
245  //Arrays for computing the track TPC using midpoints
250 
251  //Array to store output
252  double output[kMaxEvent][5][5];
253 
254  // double xyz[3];
255  // double abc[3];
256  // int chan;
257  // int cryo = fGeometry->Ncryostats();
258  // for (int c = 0; c<cryo; ++c){
259  // int tpc =fGeometry->NTPC(c);
260  // for (int t=0; t<tpc; ++t){
261  // int Nplanes=fGeometry->Nplanes(t,c);
262  // for (int p=0; p<Nplanes; ++p) {
263  // int Nwires = fGeometry->Nwires(p,t,c);
264  // for (int w=0; w<Nwires; ++w){
265  // fGeometry->WireEndPoints(c,t,p,w,xyz,abc);
266  // chan=fGeometry->PlaneWireToChannel(p,w,t,c);
267 
268 
269  // if (t % 2 != 0 && p == 2){
270  // std::cout << "FLAG " << chan << " " << c << " " << t << " " << p << " " << w << " " << xyz[0] << " " << xyz[1] << " " << xyz[2] << " " << abc[0] << " " << abc[1] << " " << abc[2] << std::endl;
271  // }
272  // //c=cryo, t=tpc, p=plane, w=wire, start (x,y,z), end (x,y,z)
273 
274  // }
275  // }
276  // }
277  // }
278 
279  //hard code (need to soft code) some detector YZ detector boundaries using geometry v5
280  double upper1 = 113.142; //Upper & lower refer to Y
281  double lower1 = -82.308;
282  double left1 = 0.29937; //Left and right refer to Z
283  double right1 = 50.1445;
284 
285  double upper3 = -1.44705;
286  double lower3 = -84.4553;
287  double left3 = 52.6724;
288  double right3 = 102.518;
289 
290  double upper5 = 113.142;
291  double lower5 = 1.46205;
292  double left5 = 52.2234;
293  double right5 = 102.068;
294 
295  double upper7 = 113.142;
296  double lower7 = -82.308;
297  double left7 = 104.147;
298  double right7 = 153.992;
299 
300 };
301 
303  : EDAnalyzer(pset)
304  , fHitsModuleLabel ( pset.get< std::string >("HitsModuleLabel"))
305  , fTrackModuleLabel ( pset.get< std::string >("TrackModuleLabel"))
306  , fClusterModuleLabel ( pset.get< std::string >("ClusterModuleLabel"))
307  , fTrkSpptAssocModuleLabel ( pset.get< std::string >("TrkSpptAssocModuleLabel"))
308  , fHitSpptAssocModuleLabel ( pset.get< std::string >("HitSpptAssocModuleLabel"))
309  , fSimulationProducerLabel ( pset.get< std::string >("SimulationLabel"))
310  , fCalorimetryModuleLabel ( pset.get< std::string >("CalorimetryModuleLabel"))
311  , fMCTruthT0ModuleLabel ( pset.get< std::string >("MCTruthT0ModuleLabel"))
312 {
313 }
314 
316 {
317  // Clean up dynamic memory and other resources here.
318 }
319 
321 {
322  // Implementation of required member function here.
323  ResetVars();
324 
326 // art::ServiceHandle<cheat::BackTrackerService> bktrk;
327 // Commenting this out since it does not appear to be used anywhere.
328 
330  art::TFileDirectory topdir = tfs->mkdir("trkgaps", "Gap histograms");
331  art::TFileDirectory gap1 = topdir.mkdir("Gap 1");
332  art::TFileDirectory gap2 = topdir.mkdir("Gap 2");
333  art::TFileDirectory gap3 = topdir.mkdir("Gap 3");
334  art::TFileDirectory gap4 = topdir.mkdir("Gap 4");
335  art::TFileDirectory gap5 = topdir.mkdir("Gap 5");
336 
337  TH1D * gapit1[kMaxEvent];
338  TH1D * gapdif1[kMaxEvent];
339 
340  TH1D * gapit2[kMaxEvent];
341  TH1D * gapdif2[kMaxEvent];
342 
343  TH1D * gapit3[kMaxEvent];
344  TH1D * gapdif3[kMaxEvent];
345 
346  TH1D * gapit4[kMaxEvent];
347  TH1D * gapdif4[kMaxEvent];
348 
349  TH1D * gapit5[kMaxEvent];
350  TH1D * gapdif5[kMaxEvent];
351 
352  run = evt.run();
353  subrun = evt.subRun();
354  event = evt.id().event();
355  art::Timestamp ts = evt.time();
356  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
357  evttime = tts.AsDouble();
358 
359  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
360  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
361  efield[0] = detProp.Efield(0);
362  efield[1] = detProp.Efield(1);
363  efield[2] = detProp.Efield(2);
364 
365  t0 = trigger_offset(clockData);
366 
367  art::Handle< std::vector<recob::Track> > trackListHandle;
368  std::vector<art::Ptr<recob::Track> > tracklist;
369  if (evt.getByLabel(fTrackModuleLabel,trackListHandle))
370  art::fill_ptr_vector(tracklist, trackListHandle);
371 
372  art::Handle< std::vector<recob::Hit> > hitListHandle;
373  std::vector<art::Ptr<recob::Hit> > hitlist;
374  if (evt.getByLabel(fHitsModuleLabel,hitListHandle))
375  art::fill_ptr_vector(hitlist, hitListHandle);
376 
377  art::Handle< std::vector<recob::Cluster> > clusterListHandle;
378  std::vector<art::Ptr<recob::Cluster> > clusterlist;
379  if (evt.getByLabel(fClusterModuleLabel,clusterListHandle))
380  art::fill_ptr_vector(clusterlist, clusterListHandle);
381 
383  evt.getByLabel(fSimulationProducerLabel, simChannelHandle);
384 
386  evt.getByLabel(fTrackModuleLabel, trackh);
387 
389  evt.getByLabel(fTrackModuleLabel, trackvh);
390 
391  // Protect against invalid art::Handle (both here and when trackvh is used below)
392  // TODO: How to do this for art::PtrVector without the uninitialised iterator?
393  std::vector< art::PtrVector<recob::Track> >::const_iterator cti;
394  if (trackvh.isValid()) cti = trackvh->begin();
395 
396  //track information
397  ntracks_reco=tracklist.size();
398 
399  recob::Track::Point_t larStart;
400  recob::Track::Point_t larEnd;
401 
402  // Get Cryostat information.....
403  int c=0;//only one cryo
404  double boundaries[6];
405  double TempBounds[6];
406  int NCryo = geom->Ncryostats();
407  for ( int b=0; b<6; b++) boundaries[b] = 0;
408  for ( int c1 = 0; c1 < NCryo; c1++ ) {
409  geom->CryostatBoundaries(TempBounds, c); // Cryostat boundaries ( neg x, pos x, neg y, pos y, neg z, pos z )
410  if ( boundaries[0] > TempBounds [0] ) boundaries[0] = TempBounds [0];
411  if ( boundaries[1] < TempBounds [1] ) boundaries[1] = TempBounds [1];
412  if ( boundaries[2] > TempBounds [2] ) boundaries[2] = TempBounds [2];
413  if ( boundaries[3] < TempBounds [3] ) boundaries[3] = TempBounds [3];
414  if ( boundaries[4] > TempBounds [4] ) boundaries[4] = TempBounds [4];
415  if ( boundaries[5] < TempBounds [5] ) boundaries[5] = TempBounds [5];
416  }
417  //std::cout << boundaries[0] << " " << boundaries[1] << " " << boundaries[2] << " " << boundaries[3] << " " <<boundaries[4] << " " << boundaries[5] << std::endl;
418 
419  // ------------------------------------
420  // NOW DO THE TRACK LIST STUFF
421  // ------------------------------------
422 
423  if ( trackListHandle.isValid() ) {
424  art::FindManyP<recob::SpacePoint> fmsp (trackListHandle, evt, fTrackModuleLabel);
425  art::FindManyP<recob::Hit> fmht (trackListHandle, evt, fTrackModuleLabel);
426  art::FindMany<anab::Calorimetry> fmcal (trackListHandle, evt, fCalorimetryModuleLabel);
427  art::FindMany<anab::T0> fmt0 (trackListHandle, evt, fMCTruthT0ModuleLabel);
428 
429  for(int i=0; i<std::min(int(tracklist.size()),kMaxTrack);++i){
430 
431 
432  //***************************************************
433  // T0 stuff - So can correct X Start/End positions
434  //***************************************************
435  double TickT0 = 0;
436  if ( fmt0.isValid() ) {
437  std::vector<const anab::T0*> T0s = fmt0.at(i);
438  for (size_t t0size =0; t0size < T0s.size(); t0size++) {
439  trkMCTruthT0[i] = T0s[t0size]->Time();
440  TickT0 = trkMCTruthT0[i] / sampling_rate(clockData);
441  } // T0 size
442  } // T0 valid
443  // Add other T0 handles...with same structure...
444  //**************
445  // END T0 stuff
446  //***************
447 
448  //******************************
449  // Hit Level Stuff - Correct X
450  //******************************
451  std::vector< art::Ptr<recob::Hit> > allHits = fmht.at(i);
452  double Hit_Size = allHits.size();
453  //*********************************
454  // End Hit Level Stuff - Correct X
455  //*********************************
456 
457  trkid[i]=i;
458  recob::Track::Point_t trackStart, trackEnd;
459  std::tie(trackStart, trackEnd) = tracklist[i]->Extent();
460  larStart = tracklist[i]->VertexDirection();
461  larEnd = tracklist[i]->EndDirection();
462 
463  trkstartx[i] = trackStart.X() - detProp.ConvertTicksToX( TickT0, allHits[Hit_Size-1]->WireID().Plane, allHits[Hit_Size-1]->WireID().TPC, allHits[Hit_Size-1]->WireID().Cryostat ); // Correct X, last entry is first 'hit'
464  trkstarty[i] = trackStart.Y();
465  trkstartz[i] = trackStart.Z();
466  trkendx[i] = trackEnd.X() - detProp.ConvertTicksToX( TickT0, allHits[0]->WireID().Plane, allHits[0]->WireID().TPC, allHits[0]->WireID().Cryostat ); // Correct X, first entry is last 'hit'
467  trkendy[i] = trackEnd.Y();
468  trkendz[i] = trackEnd.Z();
469  trkstartdcosx[i] = larStart.X();
470  trkstartdcosy[i] = larStart.Y();
471  trkstartdcosz[i] = larStart.Z();
472  trkenddcosx[i] = larEnd.X();
473  trkenddcosy[i] = larEnd.Y();
474  trkenddcosz[i] = larEnd.Z();
475  TLorentzVector v1(trackStart.X(),trackStart.Y(),trackStart.Z(),0);
476  TLorentzVector v2(trackEnd.X(),trackEnd.Y(),trackEnd.Z(),0);
477  trklen[i]=(v2-v1).Rho();
478 
479  art::Ptr<recob::Track> ptrack(trackh, i);
480  const recob::Track& track = *ptrack;
481  trklen_L[i]=track.Length();
482 
483  // Fill histograms involving reco tracks only.
484  int ntraj = track.NumberTrajectoryPoints();
485  if(ntraj > 0) {
486  auto dir = track.VertexDirection();
487  trktheta_xz[i] = std::atan2(dir.X(), dir.Z());
488  trktheta_yz[i] = std::atan2(dir.Y(), dir.Z());
489  trketa_xy[i] = std::atan2(dir.X(), dir.Y());
490  trketa_zy[i] = std::atan2(dir.Z(), dir.Y());
491  trktheta[i]=dir.Theta();
492  trkphi[i]=dir.Phi();
493  }
494 
495  double distance_squared=0;
496  TVector3 V1(trackStart.X(),trackStart.Y(),trackStart.Z());
497  TVector3 V2(trackEnd.X(),trackEnd.Y(),trackEnd.Z());
498  TVector3 vOrth=(V2-V1).Orthogonal();
499  TVector3 pointVector=V1;
500 
501  // *********************
502  // Space Point stuff:
503  // *********************
504  if(fmsp.isValid() ){
505  ntrkhits[i] = fmsp.at(i).size();
506  //std::cout << "Space Points " << ntrkhits[i] << std::endl;
507  // double distance_squared=0;
508  double distance=0;
509 
510  std::vector<art::Ptr<recob::SpacePoint> > spts = fmsp.at(i);
511  for (size_t j = 0; j<spts.size(); ++j){
512  TVector3 sptVector(spts[j]->XYZ()[0],spts[j]->XYZ()[1],spts[j]->XYZ()[2]);
513  TVector3 vToPoint=sptVector-pointVector;
514  distance=(vOrth.Dot(vToPoint))/vOrth.Mag();
515  distance_squared+=distance *distance;
516 
517  if (spts[j]->XYZ()[0] < -1000) continue;
518  else trkx[i][j] = spts[j]->XYZ()[0];
519  if (spts[j]->XYZ()[1] < -1000) continue;
520  else trky[i][j] = spts[j]->XYZ()[1];
521  if (spts[j]->XYZ()[2] < -1000) continue;
522  else trky[i][j] = spts[j]->XYZ()[2];
523 
524  std::cout << "X, Y, Z: " << spts[j]->XYZ()[0] << ", " << spts[j]->XYZ()[1] << ", " << spts[j]->XYZ()[2] << std::endl;
525 
526  }
527  distance_squared=distance_squared/spts.size();
528  trkd2[i]=distance_squared;
529 
530  }
531  // ***********************
532  // END Space Point stuff:
533  // ***********************
534 
535  // *********************
536  // Calorimetric stuff:
537  // *********************
538  if (fmcal.isValid()){
539  std::vector<const anab::Calorimetry*> calos = fmcal.at(i);
540  //std::cout<<"calo size "<<calos.size()<<std::endl;
541  for (size_t jj = 0; jj<calos.size(); ++jj){
542  trkrange[i][jj] = calos[jj]->Range();
543  trkkinE[i][jj] = (calos[jj] -> KineticEnergy());
544  //std::cout << trkkinE[i][jj] << std::endl;
545  int tt= calos[jj] -> dEdx().size();
546  for(int k = 0; k < tt; ++k) {
547  trkdedx2[i][jj][k] = (calos[jj] -> dEdx())[k];
548  trkdqdx[i][jj][k] = (calos[jj] -> dQdx())[k];
549  trkpitchHit[i][jj][k] = (calos[jj] -> TrkPitchVec())[k];
550  trkresrg[i][jj][k] = (calos[jj] -> ResidualRange())[k];
551  trkTPC[i][jj][k] = (calos[jj]->PlaneID()).TPC;
552  trkplaneid[i][jj][k] = (calos[jj]->PlaneID()).Plane;
553 
554  trkPosx[i][jj][k] = (calos[jj]->XYZ())[k].x();
555  trkPosy[i][jj][k] = (calos[jj]->XYZ())[k].y();
556  trkPosz[i][jj][k] = (calos[jj]->XYZ())[k].z();
557 
558  trkdQdxSum[i] += trkdqdx[i][jj][k];
559  trkdEdxSum[i] += trkdedx2[i][jj][k];
560  }
561  }
562  trkdQdxAverage[i] = trkdQdxSum[i] / calos.size();
563  trkdEdxAverage[i] = trkdEdxSum[i] / calos.size();
564  }
565  // ************************
566  // End Calorimetric stuff
567  // ************************
568 
569  // ---------FIND THE TRACK PITCH IN EACH PLANE------------
570  for (int j = 0; j<3; ++j){
571  try {
572  if (j==0)
573  trkpitch[i][j] = lar::util::TrackPitchInView(*(tracklist[i]), geo::kU);
574  else if (j==1)
575  trkpitch[i][j] = lar::util::TrackPitchInView(*(tracklist[i]), geo::kV);
576  else if (j==2)
577  trkpitch[i][j] = lar::util::TrackPitchInView(*(tracklist[i]), geo::kZ);
578  }
579  catch( cet::exception &e) {
580  mf::LogWarning("GapWidth")<<"caught exception "<<e<<"\n setting pitch to 0";
581  trkpitch[i][j] = 0;
582  }
583  }
584 
585  } // TRACK LOOP
586  }
587  // -----------------------------------------------
588  // NOW DO THE CLUSTER/HIT STUFF.......
589  // -----------------------------------------------
590  if ( hitListHandle.isValid() ) {
591  art::FindMany<recob::Track> fmtk (hitListHandle , evt, fTrackModuleLabel);
592  nhits = hitlist.size();
593  nhits2 = std::min(int(hitlist.size()),kMaxHits);
594  nclust = clusterlist.size();
595  for (int i = 0; i < nhits2; ++i){
596  unsigned int channel = hitlist[i]->Channel();
597  geo::WireID wireid = hitlist[i]->WireID();
598  hit_tpc[i] = wireid.TPC;
599  hit_plane[i] = wireid.Plane;
600  hit_wire[i] = wireid.Wire;
601  hit_channel[i] = channel;
602  hit_peakT[i] = hitlist[i]->PeakTime();
603  hit_charge[i] = hitlist[i]->Integral();
604  hit_ph[i] = hitlist[i]->PeakAmplitude();
605  if (fmtk.at(i).size()!=0 && fmtk.at(i)[0]->ID() <= 7){
606  std::cout << "Hit TrackID: " << fmtk.at(i)[0]->ID() << std::endl;
607  hit_trkid[i] = fmtk.at(i)[0]->ID();
608  }
609  if ( i == kMaxHits ) break;
610  }
611  }
612 
613  gapit1[event] = gap1.make<TH1D>(Form("gapit1, from event %d",event),Form("Gap Width Against Translation, Event %d",event), 101, -5.05, 5.05);
614  gapdif1[event] = gap1.make<TH1D>(Form("gapdif1 from event %d",event),Form("Gap Difference from Nominal, Event %d",event), 101, -5.05, 5.05);
615 
616  gapit2[event] = gap2.make<TH1D>(Form("gapit2, from event %d",event),Form("Gap Width Against Translation, Event %d",event), 101, -5.05, 5.05);
617  gapdif2[event] = gap2.make<TH1D>(Form("gapdif2 from event %d",event),Form("Gap Difference from Nominal, Event %d",event), 101, -5.05, 5.05);
618 
619  gapit3[event] = gap3.make<TH1D>(Form("gapit3, from event %d",event),Form("Gap Width Against Translation, Event %d",event), 101, -5.05, 5.05);
620  gapdif3[event] = gap3.make<TH1D>(Form("gapdif3 from event %d",event),Form("Gap Difference from Nominal, Event %d",event), 101, -5.05, 5.05);
621 
622  gapit4[event] = gap4.make<TH1D>(Form("gapit4, from event %d",event),Form("Gap Width Against Translation, Event %d",event), 101, -5.05, 5.05);
623  gapdif4[event] = gap4.make<TH1D>(Form("gapdif4 from event %d",event),Form("Gap Difference from Nominal, Event %d",event), 101, -5.05, 5.05);
624 
625  gapit5[event] = gap5.make<TH1D>(Form("gapit5, from event %d",event),Form("Gap Width Against Translation, Event %d",event), 101, -5.05, 5.05);
626  gapdif5[event] = gap5.make<TH1D>(Form("gapdif5 from event %d",event),Form("Gap Difference from Nominal, Event %d",event), 101, -5.05, 5.05);
627 
628  //First loop assigns tracks TPCs and limits tracks such that they don't extend beyond TPC's they occur in
629  for (int i = 0; i < ntracks_reco; i++) {
630 
631  //First, find track TPC's using mid points
632 
633  trackmidx[i] = (trkendx[i] + trkstartx[i])/2.;
634  trackmidy[i] = (trkendy[i] + trkstarty[i])/2.;
635  trackmidz[i] = (trkendz[i] + trkstartz[i])/2.;
636 
637  if (trackmidx[i] > 0.) {
638  if (trackmidz[i] > 0. && trackmidz[i] < 51.) tracktpc[i]=1;
639  if (trackmidz[i] > 51. && trackmidz[i] < 103. && trackmidy[i] < 0.) tracktpc[i]=3;
640  if (trackmidz[i] > 51. && trackmidz[i] < 103. && trackmidy[i] > 0.) tracktpc[i]=5;
641  if (trackmidz[i] > 103. && trackmidz[i] < 155.) tracktpc[i]=7;
642  }
643 
644  //Floor boundaries to coincide with wire ends
645 
646  if (tracktpc[i] == 1) {
647  if (trkendy[i] > upper1) trkendy[i] = upper1;
648  if (trkstarty[i] > upper1) trkstarty[i] = upper1;
649  if (trkendy[i] < lower1) trkendy[i] = lower1;
650  if (trkstarty[i] < lower1) trkstarty[i] = lower1;
651  if (trkendz[i] < left1 ) trkendz[i] = left1;
652  if (trkstartz[i] < left1 ) trkstartz[i] = left1;
653  if (trkendz[i] > right1) trkendz[i] = right1;
654  if (trkstartz[i] > right1) trkstartz[i] = right1;
655  }
656 
657  if (tracktpc[i] == 3) {
658  if (trkendy[i] > upper3) trkendy[i] = upper3;
659  if (trkstarty[i] > upper3) trkstarty[i] = upper3;
660  if (trkendy[i] < lower3) trkendy[i] = lower3;
661  if (trkstarty[i] < lower3) trkstarty[i] = lower3;
662  if (trkendz[i] < left3 ) trkendz[i] = left3;
663  if (trkstartz[i] < left3 ) trkstartz[i] = left3;
664  if (trkendz[i] > right3) trkendz[i] = right3;
665  if (trkstartz[i] > right3) trkstartz[i] = right3;
666  }
667 
668  if (tracktpc[i] == 5) {
669  if (trkendy[i] > upper5) trkendy[i] = upper5;
670  if (trkstarty[i] > upper5) trkstarty[i] = upper5;
671  if (trkendy[i] < lower5) trkendy[i] = lower5;
672  if (trkstarty[i] < lower5) trkstarty[i] = lower5;
673  if (trkendz[i] < left5 ) trkendz[i] = left5;
674  if (trkstartz[i] < left5 ) trkstartz[i] = left5;
675  if (trkendz[i] > right5) trkendz[i] = right5;
676  if (trkstartz[i] > right5) trkstartz[i] = right5;
677  }
678 
679  if (tracktpc[i] == 7) {
680  if (trkendy[i] > upper7) trkendy[i] = upper7;
681  if (trkstarty[i] > upper7) trkstarty[i] = upper7;
682  if (trkendy[i] < lower7) trkendy[i] = lower7;
683  if (trkstarty[i] < lower7) trkstarty[i] = lower7;
684  if (trkendz[i] < left7 ) trkendz[i] = left7;
685  if (trkstartz[i] < left7 ) trkstartz[i] = left7;
686  if (trkendz[i] > right7) trkendz[i] = right7;
687  if (trkstartz[i] > right7) trkstartz[i] = right7;
688  }
689 
690  tracklengthXZ[i] = TMath::Power((TMath::Power(std::fabs(trkendz[i] - trkstartz[i]), 2) + TMath::Power(std::fabs(trkendx[i] - trkstartx[i]), 2)), 0.5);
691  tracklengthYZ[i] = TMath::Power((TMath::Power(std::fabs(trkendz[i] - trkstartz[i]), 2) + TMath::Power(std::fabs(trkendy[i] - trkstarty[i]), 2)), 0.5);
692 
693  if (trkendz[i]>trkstartz[i]) xdiff = (trkendx[i] - trkstartx[i]);
694  else xdiff = (trkstartx[i] - trkendx[i]);
695 
696  if (trkendz[i]>trkstartz[i]) ydiff = (trkendy[i] - trkstarty[i]);
697  else ydiff = (trkstarty[i] - trkendy[i]);
698 
699  if (trkendz[i]>trkstartz[i]) zdiff = fabs(trkendz[i] - trkstartz[i]);
700  else zdiff = fabs(trkstartz[i] - trkendz[i]);
701 
702  trackgradientXZ[i] = zdiff/xdiff;
703  trackgradientYZ[i] = zdiff/ydiff;
704  }
705 
706  for (int translationsteps = 0; translationsteps < 101; translationsteps ++){
707  translation = (double(translationsteps)-50)/10.;
708 
709  //Start of loop for track comparison
710  for (int i = 0; i < ntracks_reco; i++){
711 
712  //need comparison track
713  for (int k = 0; k < ntracks_reco; k++){
714 
715  possiblex.clear();
716  possibley.clear();
717  possiblez.clear();
718 
719  track1x = 0;
720  track1y = 0;
721  track1z = 0;
722  track2x = 0;
723  track2y = 0;
724  track2z = 0;
725  track1grad = 0;
726  track2grad = 0;
727  posminx = 0;
728  posminy = 0;
729  posminz = 0;
730 
731  //force tracks to be adjacent in tpcs with similar gradients
732  if ((tracktpc[i] != tracktpc[k]) && (abs(tracktpc[i] - tracktpc[k]) <= 4) && (fabs(atan(trackgradientXZ[k]))*TMath::Pi()/180 - (atan(trackgradientXZ[i]))*TMath::Pi()/180 < 0.1) && (tracktpc[i] < tracktpc[k])) {
733 
734  possiblex.push_back(fabs(trkendx[i] - trkendx[k]));
735  possiblex.push_back(fabs(trkendx[i] - trkstartx[k]));
736  possiblex.push_back(fabs(trkstartx[i] - trkendx[k]));
737  possiblex.push_back(fabs(trkstartx[i] - trkstartx[k]));
738 
739  possibley.push_back(fabs(trkendy[i] - trkendy[k]));
740  possibley.push_back(fabs(trkendy[i] - trkstarty[k]));
741  possibley.push_back(fabs(trkstarty[i] - trkendy[k]));
742  possibley.push_back(fabs(trkstarty[i] - trkstarty[k]));
743 
744  possiblez.push_back(fabs(trkendz[i] - trkendz[k]));
745  possiblez.push_back(fabs(trkendz[i] - trkstartz[k]));
746  possiblez.push_back(fabs(trkstartz[i] - trkendz[k]));
747  possiblez.push_back(fabs(trkstartz[i] - trkstartz[k]));
748 
749  minx = * std::min_element(possiblex.begin(), possiblex.end());
750  miny = * std::min_element(possibley.begin(), possibley.end());
751  minz = * std::min_element(possiblez.begin(), possiblez.end());
752 
753  for (size_t g = 0; g < possiblex.size(); g++) {
754  if (possiblex[g] == minx) {
755  posminx = g;
756  }
757  }
758  for (size_t g = 0; g < possibley.size(); g++) {
759  if (possibley[g] == miny) {
760  posminy = g;
761  }
762  }
763  for (size_t g = 0; g < possiblez.size(); g++) {
764 
765  if (possiblez[g] == minz) {
766  posminz = g;
767  }
768  }
769 
770 
771  //GAP1
772  if ((tracktpc[i] == 1 && tracktpc[k] == 5) || (tracktpc[i] == 5 && tracktpc[k] == 1)) {
773  if (minz <= 6.) {
774  //std::cout << "From the comparison between tracks with IDs: " << trkid[i] << " and " << trkid[k];
775  //std::cout << ", the size of Gap 1 is: " << minz << " at position " << posminz << std::endl;
776 
777 
778  if (posminz <= 1 && posminz % 2 ==1){
779  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
780 
781  if (trkendz[i] < trkstartz[k]){
782  track1x = trkendx[i];
783  track1y = trkendy[i];
784  track1z = trkendz[i];
785  track2x = trkstartx[k] + translation;
786  track2y = trkstarty[k];
787  track2z = trkstartz[k];
788  track1grad = trackgradientXZ[i];
789  }
790  else {
791  track2x = trkendx[i] + translation;
792  track2y = trkendy[i];
793  track2z = trkendz[i];
794  track1x = trkstartx[k];
795  track1y = trkstarty[k];
796  track1z = trkstartz[k];
797  track1grad = trackgradientXZ[k];
798  }
799  }
800 
801  if (posminz <= 1 && posminz % 2 !=1){
802  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
803 
804  if (trkendz[i] < trkendz[k]){
805  track1x = trkendx[i];
806  track1y = trkendy[i];
807  track1z = trkendz[i];
808  track2x = trkendx[k] + translation;
809  track2y = trkendy[k];
810  track2z = trkendz[k];
811  track1grad = trackgradientXZ[i];
812  }
813  else {
814  track2x = trkendx[i] + translation;
815  track2y = trkendy[i];
816  track2z = trkendz[i];
817  track1x = trkendx[k];
818  track1y = trkendy[k];
819  track1z = trkendz[k];
820  track1grad = trackgradientXZ[k];
821  }
822  }
823 
824  if (posminz > 1 && posminz % 2 ==1){
825  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
826 
827  if (trkstartz[i] < trkstartz[k]){
828  track1x = trkstartx[i];
829  track1y = trkstarty[i];
830  track1z = trkstartz[i];
831  track2x = trkstartx[k] + translation;
832  track2y = trkstarty[k];
833  track2z = trkstartz[k];
834  track1grad = trackgradientXZ[i];
835  }
836  else {
837  track2x = trkstartx[i] + translation;
838  track2y = trkstarty[i];
839  track2z = trkstartz[i];
840  track1x = trkstartx[k];
841  track1y = trkstarty[k];
842  track1z = trkstartz[k];
843  track1grad = trackgradientXZ[k];
844  }
845  }
846  if (posminz > 1 && posminz % 2 !=1){
847  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
848 
849  if (trkstartz[i] < trkendz[k]){
850  track1x = trkstartx[i];
851  track1y = trkstarty[i];
852  track1z = trkstartz[i];
853  track2x = trkendx[k] + translation;
854  track2y = trkendy[k];
855  track2z = trkendz[k];
856  track1grad = trackgradientXZ[i];
857  }
858  else {
859  track2x = trkstartx[i] + translation;
860  track2y = trkstarty[i];
861  track2z = trkstartz[i];
862  track1x = trkendx[k];
863  track1y = trkendy[k];
864  track1z = trkendz[k];
865  track1grad = trackgradientXZ[k];
866  }
867  }
868 
869  double error = 0.01;
870  TMinuit* mingen = new TMinuit(1);
871  mingen->SetFCN(minuitFunctionx);
872  mingen->SetPrintLevel(-1);
873  mingen->DefineParameter(0,"gap extension", -20, 0.001, -100, 100);
874  mingen->Migrad();
875  mingen->GetParameter(0, z, error);
876 
877  output[event][0][0] = double(trkid[i]);
878  output[event][0][1] = double(trkid[k]);
879  output[event][0][2] = double(tracktpc[i]);
880  output[event][0][3] = double(tracktpc[k]);
881  output[event][0][4] = (minz + z);
882 
883  }
884  } // End Gap 1
885 
886  //GAP2
887  if ((tracktpc[i] == 5 && tracktpc[k] == 7) || (tracktpc[i] == 7 && tracktpc[k] == 5)) {
888  if (minz <= 6.) {
889 
890  if (posminz <= 1 && posminz % 2 ==1){
891  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
892 
893  if (trkendz[i] < trkstartz[k]){
894  track1x = trkendx[i];
895  track1y = trkendy[i];
896  track1z = trkendz[i];
897  track2x = trkstartx[k] + translation;
898  track2y = trkstarty[k];
899  track2z = trkstartz[k];
900  track1grad = trackgradientXZ[i];
901  }
902  else {
903  track2x = trkendx[i] + translation;
904  track2y = trkendy[i];
905  track2z = trkendz[i];
906  track1x = trkstartx[k];
907  track1y = trkstarty[k];
908  track1z = trkstartz[k];
909  track1grad = trackgradientXZ[k];
910  }
911  }
912 
913  if (posminz <= 1 && posminz % 2 !=1){
914  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
915 
916  if (trkendz[i] < trkendz[k]){
917  track1x = trkendx[i];
918  track1y = trkendy[i];
919  track1z = trkendz[i];
920  track2x = trkendx[k] + translation;
921  track2y = trkendy[k];
922  track2z = trkendz[k];
923  track1grad = trackgradientXZ[i];
924  }
925  else {
926  track2x = trkendx[i] + translation;
927  track2y = trkendy[i];
928  track2z = trkendz[i];
929  track1x = trkendx[k];
930  track1y = trkendy[k];
931  track1z = trkendz[k];
932  track1grad = trackgradientXZ[k];
933  }
934  }
935 
936  if (posminz > 1 && posminz % 2 ==1){
937  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
938 
939  if (trkstartz[i] < trkstartz[k]){
940  track1x = trkstartx[i];
941  track1y = trkstarty[i];
942  track1z = trkstartz[i];
943  track2x = trkstartx[k] + translation;
944  track2y = trkstarty[k];
945  track2z = trkstartz[k];
946  track1grad = trackgradientXZ[i];
947  }
948  else {
949  track2x = trkstartx[i] + translation;
950  track2y = trkstarty[i];
951  track2z = trkstartz[i];
952  track1x = trkstartx[k];
953  track1y = trkstarty[k];
954  track1z = trkstartz[k];
955  track1grad = trackgradientXZ[k];
956  }
957  }
958  if (posminz > 1 && posminz % 2 !=1){
959  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
960 
961  if (trkstartz[i] < trkendz[k]){
962  track1x = trkstartx[i];
963  track1y = trkstarty[i];
964  track1z = trkstartz[i];
965  track2x = trkendx[k] + translation;
966  track2y = trkendy[k];
967  track2z = trkendz[k];
968  track1grad = trackgradientXZ[i];
969  }
970  else {
971  track2x = trkstartx[i] + translation;
972  track2y = trkstarty[i];
973  track2z = trkstartz[i];
974  track1x = trkendx[k];
975  track1y = trkendy[k];
976  track1z = trkendz[k];
977  track1grad = trackgradientXZ[k];
978  }
979  }
980 
981  double error = 0.01;
982  TMinuit* mingen = new TMinuit(1);
983  mingen->SetFCN(minuitFunctionx);
984  mingen->SetPrintLevel(-1);
985  mingen->DefineParameter(0,"gap extension", -20, 0.001, -100, 100);
986  mingen->Migrad();
987  mingen->GetParameter(0, z, error);
988 
989  output[event][1][0] = double(trkid[i]);
990  output[event][1][1] = double(trkid[k]);
991  output[event][1][2] = double(tracktpc[i]);
992  output[event][1][3] = double(tracktpc[k]);
993  output[event][1][4] = (minz + z);
994 
995  }
996  } // End Gap 2
997 
998  //GAP3
999  if ((tracktpc[i] == 1 && tracktpc[k] == 3) || (tracktpc[i] == 3 && tracktpc[k] == 1)) {
1000  if (minz <= 6.) {
1001 
1002  if (posminz <= 1 && posminz % 2 ==1){
1003  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
1004 
1005  if (trkendz[i] < trkstartz[k]){
1006  track1x = trkendx[i];
1007  track1y = trkendy[i];
1008  track1z = trkendz[i];
1009  track2x = trkstartx[k] + translation;
1010  track2y = trkstarty[k];
1011  track2z = trkstartz[k];
1012  track1grad = trackgradientXZ[i];
1013  }
1014  else {
1015  track2x = trkendx[i] + translation;
1016  track2y = trkendy[i];
1017  track2z = trkendz[i];
1018  track1x = trkstartx[k];
1019  track1y = trkstarty[k];
1020  track1z = trkstartz[k];
1021  track1grad = trackgradientXZ[k];
1022  }
1023  }
1024 
1025  if (posminz <= 1 && posminz % 2 !=1){
1026  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
1027 
1028  if (trkendz[i] < trkendz[k]){
1029  track1x = trkendx[i];
1030  track1y = trkendy[i];
1031  track1z = trkendz[i];
1032  track2x = trkendx[k] + translation;
1033  track2y = trkendy[k];
1034  track2z = trkendz[k];
1035  track1grad = trackgradientXZ[i];
1036  }
1037  else {
1038  track2x = trkendx[i] + translation;
1039  track2y = trkendy[i];
1040  track2z = trkendz[i];
1041  track1x = trkendx[k];
1042  track1y = trkendy[k];
1043  track1z = trkendz[k];
1044  track1grad = trackgradientXZ[k];
1045  }
1046  }
1047 
1048  if (posminz > 1 && posminz % 2 ==1){
1049  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
1050 
1051  if (trkstartz[i] < trkstartz[k]){
1052  track1x = trkstartx[i];
1053  track1y = trkstarty[i];
1054  track1z = trkstartz[i];
1055  track2x = trkstartx[k] + translation;
1056  track2y = trkstarty[k];
1057  track2z = trkstartz[k];
1058  track1grad = trackgradientXZ[i];
1059  }
1060  else {
1061  track2x = trkstartx[i] + translation;
1062  track2y = trkstarty[i];
1063  track2z = trkstartz[i];
1064  track1x = trkstartx[k];
1065  track1y = trkstarty[k];
1066  track1z = trkstartz[k];
1067  track1grad = trackgradientXZ[k];
1068  }
1069  }
1070  if (posminz > 1 && posminz % 2 !=1){
1071  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
1072 
1073  if (trkstartz[i] < trkendz[k]){
1074  track1x = trkstartx[i];
1075  track1y = trkstarty[i];
1076  track1z = trkstartz[i];
1077  track2x = trkendx[k] + translation;
1078  track2y = trkendy[k];
1079  track2z = trkendz[k];
1080  track1grad = trackgradientXZ[i];
1081  }
1082  else {
1083  track2x = trkstartx[i] + translation;
1084  track2y = trkstarty[i];
1085  track2z = trkstartz[i];
1086  track1x = trkendx[k];
1087  track1y = trkendy[k];
1088  track1z = trkendz[k];
1089  track1grad = trackgradientXZ[k];
1090  }
1091  }
1092 
1093  double error = 0.01;
1094  TMinuit* mingen = new TMinuit(1);
1095  mingen->SetFCN(minuitFunctionx);
1096  mingen->SetPrintLevel(-1);
1097  mingen->DefineParameter(0,"gap extension", -20, 0.001, -100, 100);
1098  mingen->Migrad();
1099  mingen->GetParameter(0, z, error);
1100 
1101  output[event][2][0] = double(trkid[i]);
1102  output[event][2][1] = double(trkid[k]);
1103  output[event][2][2] = double(tracktpc[i]);
1104  output[event][2][3] = double(tracktpc[k]);
1105  output[event][2][4] = (minz + z);
1106 
1107  }
1108  } // End Gap 3
1109 
1110  //GAP4
1111  if ((tracktpc[i] == 3 && tracktpc[k] == 7) || (tracktpc[i] == 7 && tracktpc[k] == 3)) {
1112  if (minz <= 6.) {
1113 
1114  if (posminz <= 1 && posminz % 2 ==1){
1115  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
1116 
1117  if (trkendz[i] < trkstartz[k]){
1118  track1x = trkendx[i];
1119  track1y = trkendy[i];
1120  track1z = trkendz[i];
1121  track2x = trkstartx[k] + translation;
1122  track2y = trkstarty[k];
1123  track2z = trkstartz[k];
1124  track1grad = trackgradientXZ[i];
1125  }
1126  else {
1127  track2x = trkendx[i] + translation;
1128  track2y = trkendy[i];
1129  track2z = trkendz[i];
1130  track1x = trkstartx[k];
1131  track1y = trkstarty[k];
1132  track1z = trkstartz[k];
1133  track1grad = trackgradientXZ[k];
1134  }
1135  }
1136 
1137  if (posminz <= 1 && posminz % 2 !=1){
1138  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
1139 
1140  if (trkendz[i] < trkendz[k]){
1141  track1x = trkendx[i];
1142  track1y = trkendy[i];
1143  track1z = trkendz[i];
1144  track2x = trkendx[k] + translation;
1145  track2y = trkendy[k];
1146  track2z = trkendz[k];
1147  track1grad = trackgradientXZ[i];
1148  }
1149  else {
1150  track2x = trkendx[i] + translation;
1151  track2y = trkendy[i];
1152  track2z = trkendz[i];
1153  track1x = trkendx[k];
1154  track1y = trkendy[k];
1155  track1z = trkendz[k];
1156  track1grad = trackgradientXZ[k];
1157  }
1158  }
1159 
1160  if (posminz > 1 && posminz % 2 ==1){
1161  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
1162 
1163  if (trkstartz[i] < trkstartz[k]){
1164  track1x = trkstartx[i];
1165  track1y = trkstarty[i];
1166  track1z = trkstartz[i];
1167  track2x = trkstartx[k] + translation;
1168  track2y = trkstarty[k];
1169  track2z = trkstartz[k];
1170  track1grad = trackgradientXZ[i];
1171  }
1172  else {
1173  track2x = trkstartx[i] + translation;
1174  track2y = trkstarty[i];
1175  track2z = trkstartz[i];
1176  track1x = trkstartx[k];
1177  track1y = trkstarty[k];
1178  track1z = trkstartz[k];
1179  track1grad = trackgradientXZ[k];
1180  }
1181  }
1182  if (posminz > 1 && posminz % 2 !=1){
1183  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
1184 
1185  if (trkstartz[i] < trkendz[k]){
1186  track1x = trkstartx[i];
1187  track1y = trkstarty[i];
1188  track1z = trkstartz[i];
1189  track2x = trkendx[k] + translation;
1190  track2y = trkendy[k];
1191  track2z = trkendz[k];
1192  track1grad = trackgradientXZ[i];
1193  }
1194  else {
1195  track2x = trkstartx[i] + translation;
1196  track2y = trkstarty[i];
1197  track2z = trkstartz[i];
1198  track1x = trkendx[k];
1199  track1y = trkendy[k];
1200  track1z = trkendz[k];
1201  track1grad = trackgradientXZ[k];
1202  }
1203  }
1204 
1205  double error = 0.01;
1206  TMinuit* mingen = new TMinuit(1);
1207  mingen->SetFCN(minuitFunctionx);
1208  mingen->SetPrintLevel(-1);
1209  mingen->DefineParameter(0,"gap extension", -20, 0.001, -100, 100);
1210  mingen->Migrad();
1211  mingen->GetParameter(0, z, error);
1212 
1213  output[event][3][0] = double(trkid[i]);
1214  output[event][3][1] = double(trkid[k]);
1215  output[event][3][2] = double(tracktpc[i]);
1216  output[event][3][3] = double(tracktpc[k]);
1217  output[event][3][4] = (minz + z);
1218 
1219  }
1220  } // End Gap 4
1221 
1222  //GAP 5
1223  if ((tracktpc[i] == 3 && tracktpc[k] == 5) || (tracktpc[i] == 5 && tracktpc[k] == 3)) {
1224  if (miny <= 6.) {
1225 
1226  if (posminz <= 1 && posminz % 2 ==1){
1227  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
1228 
1229  if (trkendz[i] < trkstartz[k]){
1230  track1x = trkendx[i];
1231  track1y = trkendy[i];
1232  track1z = trkendz[i];
1233  track2x = trkstartx[k];
1234  track2y = trkstarty[k];
1235  track2z = trkstartz[k] + translation;
1236  track1grad = trackgradientYZ[i];
1237  track2grad = trackgradientYZ[k];
1238  }
1239  else {
1240  track2x = trkendx[i];
1241  track2y = trkendy[i];
1242  track2z = trkendz[i] + translation;
1243  track1x = trkstartx[k];
1244  track1y = trkstarty[k];
1245  track1z = trkstartz[k];
1246  track1grad = trackgradientYZ[k];
1247  track2grad = trackgradientYZ[i];
1248  }
1249  }
1250 
1251  if (posminz <= 1 && posminz % 2 !=1){
1252  //std::cout << "Z minimised matching end of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
1253 
1254  if (trkendz[i] < trkendz[k]){
1255  track1x = trkendx[i];
1256  track1y = trkendy[i];
1257  track1z = trkendz[i];
1258  track2x = trkendx[k];
1259  track2y = trkendy[k];
1260  track2z = trkendz[k] + translation;
1261  track1grad = trackgradientYZ[i];
1262  track2grad = trackgradientYZ[k];
1263  }
1264  else {
1265  track2x = trkendx[i];
1266  track2y = trkendy[i];
1267  track2z = trkendz[i] + translation;
1268  track1x = trkendx[k];
1269  track1y = trkendy[k];
1270  track1z = trkendz[k];
1271  track1grad = trackgradientYZ[k];
1272  track2grad = trackgradientYZ[i];
1273  }
1274  }
1275 
1276  if (posminz > 1 && posminz % 2 ==1){
1277  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the start of track: " << trkid[k] << std::endl;
1278 
1279  if (trkstartz[i] < trkstartz[k]){
1280  track1x = trkstartx[i];
1281  track1y = trkstarty[i];
1282  track1z = trkstartz[i];
1283  track2x = trkstartx[k];
1284  track2y = trkstarty[k];
1285  track2z = trkstartz[k] + translation;
1286  track1grad = trackgradientYZ[i];
1287  track2grad = trackgradientYZ[k];
1288  }
1289  else {
1290  track2x = trkstartx[i];
1291  track2y = trkstarty[i];
1292  track2z = trkstartz[i] + translation;
1293  track1x = trkstartx[k];
1294  track1y = trkstarty[k];
1295  track1z = trkstartz[k];
1296  track1grad = trackgradientYZ[k];
1297  track2grad = trackgradientYZ[i];
1298  }
1299  }
1300  if (posminz > 1 && posminz % 2 !=1){
1301  //std::cout << "Z minimised matching start of track: " << trkid[i] << " and the end of track: " << trkid[k] << std::endl;
1302 
1303  if (trkstartz[i] < trkendz[k]){
1304  track1x = trkstartx[i];
1305  track1y = trkstarty[i];
1306  track1z = trkstartz[i];
1307  track2x = trkendx[k];
1308  track2y = trkendy[k];
1309  track2z = trkendz[k] + translation;
1310  track1grad = trackgradientYZ[i];
1311  track2grad = trackgradientYZ[k];
1312  }
1313  else {
1314  track2x = trkstartx[i];
1315  track2y = trkstarty[i];
1316  track2z = trkstartz[i] + translation;
1317  track1x = trkendx[k];
1318  track1y = trkendy[k];
1319  track1z = trkendz[k];
1320  track1grad = trackgradientYZ[k];
1321  track2grad = trackgradientYZ[i];
1322  }
1323  }
1324 
1325  double error = 0.01;
1326  TMinuit* min5 = new TMinuit(1);
1327  min5->SetFCN(minuitFunctionz);
1328  min5->SetPrintLevel(-1);
1329  min5->DefineParameter(0,"gap extension", -20, 0.001, -100, 100);
1330  min5->Migrad();
1331  min5->GetParameter(0, y, error);
1332 
1333  output[event][4][0] = double(trkid[i]);
1334  output[event][4][1] = double(trkid[k]);
1335  output[event][4][2] = double(tracktpc[i]);
1336  output[event][4][3] = double(tracktpc[k]);
1337  output[event][4][4] = (miny + y);
1338 
1339  }
1340  }//End Gap 5
1341 
1342 
1343 
1344 
1345  }//End of check that TPCs are adjacent and that track ID's are not identical
1346  }//End of loop over tracks[k]
1347  } //End of loop over tracks[i]
1348 
1349  gapit1[event]->SetBinContent(translationsteps+1, output[event][0][4]);
1350  gapdif1[event]->SetBinContent(translationsteps+1, fabs(output[event][0][4] - 2.0789));
1351 
1352  gapit2[event]->SetBinContent(translationsteps+1, output[event][1][4]);
1353  gapdif2[event]->SetBinContent(translationsteps+1, fabs(output[event][1][4] - 2.079));
1354 
1355  gapit3[event]->SetBinContent(translationsteps+1, output[event][2][4]);
1356  gapdif3[event]->SetBinContent(translationsteps+1, fabs(output[event][2][4] - 2.5279));
1357 
1358  gapit4[event]->SetBinContent(translationsteps+1, output[event][3][4]);
1359  gapdif4[event]->SetBinContent(translationsteps+1, fabs(output[event][3][4] - 1.629));
1360 
1361  gapit5[event]->SetBinContent(translationsteps+1, output[event][4][4]);
1362  gapdif5[event]->SetBinContent(translationsteps+1, fabs(output[event][4][4] - 2.9091));
1363 
1364  }
1365 
1366  // -------------------------------------------
1367  // FINALLLY Fill Tree:
1368  // -------------------------------------------
1369  fTree->Fill();
1370 
1371 }
1372 
1373 
1375 {
1376  // Implementation of optional member function here.
1378  fTree = tfs2->make<TTree>("gapwidth","analysis tree");
1379 
1380  fTree->Branch("run",&run,"run/I");
1381  fTree->Branch("subrun",&subrun,"subrun/I");
1382  fTree->Branch("event",&event,"event/I");
1383  fTree->Branch("evttime",&evttime,"evttime/D");
1384  fTree->Branch("efield",efield,"efield[3]/D");
1385  fTree->Branch("t0",&t0,"t0/I");
1386 
1387  fTree->Branch("ntracks_reco",&ntracks_reco,"ntracks_reco/I");
1388  fTree->Branch("ntrkhits",ntrkhits,"ntrkhits[ntracks_reco]/I");
1389  fTree->Branch("trkid",trkid,"trkid[ntracks_reco]/I");
1390  fTree->Branch("trkstartx",trkstartx,"trkstartx[ntracks_reco]/D");
1391  fTree->Branch("trkstarty",trkstarty,"trkstarty[ntracks_reco]/D");
1392  fTree->Branch("trkstartz",trkstartz,"trkstartz[ntracks_reco]/D");
1393  fTree->Branch("trkendx",trkendx,"trkendx[ntracks_reco]/D");
1394  fTree->Branch("trkendy",trkendy,"trkendy[ntracks_reco]/D");
1395  fTree->Branch("trkendz",trkendz,"trkendz[ntracks_reco]/D");
1396  fTree->Branch("trkstartdcosx",trkstartdcosx,"trkstartdcosx[ntracks_reco]/D");
1397  fTree->Branch("trkstartdcosy",trkstartdcosy,"trkstartdcosy[ntracks_reco]/D");
1398  fTree->Branch("trkstartdcosz",trkstartdcosz,"trkstartdcosz[ntracks_reco]/D");
1399  fTree->Branch("trkenddcosx",trkenddcosx,"trkenddcosx[ntracks_reco]/D");
1400  fTree->Branch("trkenddcosy",trkenddcosy,"trkenddcosy[ntracks_reco]/D");
1401  fTree->Branch("trkenddcosz",trkenddcosz,"trkenddcosz[ntracks_reco]/D");
1402 
1403  //Array branches start
1404  fTree->Branch("trkx",trkx,"trkx[ntracks_reco][1000]/D");
1405  fTree->Branch("trky",trky,"trky[ntracks_reco][1000]/D");
1406  fTree->Branch("trkz",trkz,"trkz[ntracks_reco][1000]/D");
1407 
1408  //Array branches end
1409 
1410  fTree->Branch("trktheta_xz",trktheta_xz,"trktheta_xz[ntracks_reco]/D");
1411  fTree->Branch("trktheta_yz",trktheta_yz,"trktheta_yz[ntracks_reco]/D");
1412  fTree->Branch("trketa_xy",trketa_xy,"trketa_xy[ntracks_reco]/D");
1413  fTree->Branch("trketa_zy",trketa_zy,"trketa_zy[ntracks_reco]/D");
1414  fTree->Branch("trktheta",trktheta,"trktheta[ntracks_reco]/D");
1415  fTree->Branch("trkphi",trkphi,"trkphi[ntracks_reco]/D");
1416  fTree->Branch("trkd2",trkd2,"trkd2[ntracks_reco]/D");
1417 
1418  //Array Branches
1419  fTree->Branch("trkdedx2",trkdedx2,"trkdedx2[ntracks_reco][3][1000]/D");
1420  fTree->Branch("trkdqdx",trkdqdx,"trkdqdx[ntracks_reco][3][1000]/D");
1421 
1422  fTree->Branch("trkpitch",trkpitch,"trkpitch[ntracks_reco][3]/D");
1423 
1424  //Array
1425  fTree->Branch("trkpitchHit",trkpitchHit,"trkpitchHit[ntracks_reco][3][1000]/D");
1426  fTree->Branch("trkkinE",trkkinE,"trkkinE[ntracks_reco][3]/D");
1427  fTree->Branch("trkrange",trkrange,"trkrange[ntracks_reco][3]/D");
1428 
1429  //Array branches
1430  fTree->Branch("trkTPC",trkTPC,"trkTPC[ntracks_reco][3][1000]/D");
1431  fTree->Branch("trkplaneid",trkplaneid,"trkplaneid[ntracks_reco][3][1000]/D");
1432  fTree->Branch("trkresrg",trkresrg,"trkresrg[ntracks_reco][3][1000]/D");
1433  fTree->Branch("trkPosx",trkPosx,"trkPosx[ntracks_reco][3][1000]/D");
1434  fTree->Branch("trkPosy",trkPosy,"trkPosy[ntracks_reco][3][1000]/D");
1435  fTree->Branch("trkPosz",trkPosz,"trkPosz[ntracks_reco][3][1000]/D");
1436  //Array branches end
1437 
1438  fTree->Branch("trklen",trklen,"trklen[ntracks_reco]/D");
1439  fTree->Branch("trklen_L",trklen_L,"trklen_L[ntracks_reco]/D");
1440  fTree->Branch("trkdQdxSum",trkdQdxSum,"trkdQdxSum[ntracks_reco]/D");
1441  fTree->Branch("trkdQdxAverage",trkdQdxAverage,"trkdQdxAverage[ntracks_reco]/D");
1442  fTree->Branch("trkdEdxSum",trkdEdxSum,"trkdEdxSum[ntracks_reco]/D");
1443  fTree->Branch("trkdEdxAverage",trkdEdxAverage,"trkdEdxAverage[ntracks_reco]/D");
1444 
1445  fTree->Branch("trkMCTruthT0",trkMCTruthT0,"trkMCTruthT0[ntracks_reco]/D");
1446 
1447  fTree->Branch("nhits",&nhits,"nhits/I");
1448  fTree->Branch("nhits2",&nhits2,"nhits2/I");
1449  fTree->Branch("nclust",&nclust,"nclust/I");
1450  fTree->Branch("hit_plane",hit_plane,"hit_plane[nhits2]/I");
1451  fTree->Branch("hit_tpc",hit_tpc,"hit_tpc[nhits2]/I");
1452  fTree->Branch("hit_wire",hit_wire,"hit_wire[nhits2]/I");
1453  fTree->Branch("hit_channel",hit_channel,"hit_channel[nhits2]/I");
1454  fTree->Branch("hit_peakT",hit_peakT,"hit_peakT[nhits2]/D");
1455  fTree->Branch("hit_charge",hit_charge,"hit_charge[nhits2]/D");
1456  fTree->Branch("hit_ph",hit_ph,"hit_ph[nhits2]/D");
1457  fTree->Branch("hit_trkid",hit_trkid,"hit_trkid[nhits2]/I");
1458 
1459 }
1460 
1461 //void GapWidth::GapWidth::reconfigure(fhicl::ParameterSet const & p)
1462 //{
1463 // // Implementation of optional member function here.
1464 //}
1466 
1467  run = -99999;
1468  subrun = -99999;
1469  event = -99999;
1470  evttime = -99999;
1471  for (int i = 0; i<3; ++i){
1472  efield[i] = -99999;
1473  }
1474  t0 = -99999;
1475  ntracks_reco = -99999;
1476 
1477  for (int i = 0; i < kMaxTrack; ++i){
1478  trkstartx[i] = -99999;
1479  trkstarty[i] = -99999;
1480  trkstartz[i] = -99999;
1481  trkendx[i] = -99999;
1482  trkendy[i] = -99999;
1483  trkendz[i] = -99999;
1484 
1485  trkd2[i] = -99999;
1486  trktheta_xz[i] = -99999;
1487  trktheta_yz[i] = -99999;
1488  trketa_xy[i] = -99999;
1489  trketa_zy[i] = -99999;
1490  trktheta[i] = -99999;
1491  trkphi[i] = -99999;
1492 
1493  trkdQdxSum[i] = 0;
1494  trkdEdxSum[i] = 0;
1495  for(int ii=0;ii<3;ii++)
1496  {
1497  trkkinE[i][ii] = -99999;
1498  trkrange[i][ii] = -99999;
1499  for(int k=0;k<1000;k++)
1500  {
1501  trkdedx2[i][ii][k] = -99999;
1502  trkdqdx[i][ii][k] = -99999;
1503  trkpitchHit[i][ii][k] = -99999;
1504  trkTPC[i][ii][k] = -99999;
1505  trkplaneid[i][ii][k] = -99999;
1506  trkresrg[i][ii][k] = -99999;
1507  trkPosx[i][ii][k] = -99999;
1508  trkPosy[i][ii][k] = -99999;
1509  trkPosz[i][ii][k] = -99999;
1510  }
1511  }
1512  trkstartdcosx[i] = -99999;
1513  trkstartdcosy[i] = -99999;
1514  trkstartdcosz[i] = -99999;
1515  trklen[i] = -99999;
1516  trklen_L[i] = -99999;
1517  trkid[i] = -99999;
1518 
1519  trkenddcosx[i] = -99999;
1520  trkenddcosy[i] = -99999;
1521  trkenddcosz[i] = -99999;
1522  ntrkhits[i] = -99999;
1523  for (int j = 0; j<kMaxTrackHits; ++j){
1524  trkx[i][j] = -99999;
1525  trky[i][j] = -99999;
1526  trkz[i][j] = -99999;
1527  }
1528  for (int j = 0; j<3; ++j){
1529  trkpitch[i][j] = -99999;
1530  }
1531  }
1532  nhits = -99999;
1533  nhits2= 0;
1534  for (int i = 0; i<kMaxHits; ++i){
1535  hit_plane[i] = -99999;
1536  hit_tpc[i] = -99999;
1537  hit_wire[i] = -99999;
1538  hit_channel[i] = -99999;
1539  hit_peakT[i] = -99999;
1540  hit_charge[i] = -99999;
1541  hit_ph[i] = -99999;
1542  hit_trkid[i] = -99999;
1543  }
1544 }
1545 
1547 {
1548 }
1549 
code to link reconstructed objects back to the MC truth information
run
Definition: dbjson.py:15
double trky[kMaxTrack][kMaxTrackHits]
double tracklengthXZ[kMaxTrack]
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
double trklen_L[kMaxTrack]
std::string fCalorimetryModuleLabel
const int kMaxHits
static QCString result
double track2grad
Encapsulate the construction of a single cyostat.
const int kMaxTrackHits
double miny
double trkdQdxSum[kMaxTrack]
double trketa_zy[kMaxTrack]
double trkPosy[kMaxTrack][3][1000]
std::string string
Definition: nybbler.cc:12
Planes which measure V.
Definition: geo_types.h:130
double trkenddcosy[kMaxTrack]
double trackmidy[kMaxTrack]
GapWidth(fhicl::ParameterSet const &p)
double trkendx[kMaxTrack]
int tracktpc[kMaxTrack]
double trkPosz[kMaxTrack][3][1000]
Declaration of signal hit object.
double trackmidz[kMaxTrack]
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
double trkdEdxAverage[kMaxTrack]
std::string fClusterModuleLabel
int posminz
int hit_tpc[kMaxHits]
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
std::string fSimulationProducerLabel
std::string fTrackModuleLabel
Vector_t VertexDirection() const
Definition: Track.h:132
error
Definition: include.cc:26
STL namespace.
std::vector< double > possiblez
Planes which measure Z direction.
Definition: geo_types.h:132
double xextrap(double z)
double track2x
int trkid[kMaxTrack]
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
intermediate_table::const_iterator const_iterator
double hit_charge[kMaxHits]
static const double g
Definition: Units.h:144
double track1grad
double track2z
double track1z
string dir
double trkTPC[kMaxTrack][3][1000]
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:27
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Particle class.
double trktheta_yz[kMaxTrack]
double tracklengthYZ[kMaxTrack]
void CryostatBoundaries(double *boundaries, geo::CryostatID const &cid) const
Returns the boundaries of the specified cryostat.
std::vector< double > possibley
double trkstartx[kMaxTrack]
std::string fHitSpptAssocModuleLabel
double trkdedx2[kMaxTrack][3][1000]
double minz
double trkdqdx[kMaxTrack][3][1000]
int ntrkhits[kMaxTrack]
bool isValid() const
Definition: Handle.h:183
double trackmidx[kMaxTrack]
double y
Definition: type_traits.h:61
T abs(T value)
Planes which measure U.
Definition: geo_types.h:129
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
double dEdx(float dqdx, float Efield)
Definition: doAna.cpp:21
void beginJob() override
double trkkinE[kMaxTrack][3]
double trkphi[kMaxTrack]
const double e
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
int posminy
Timestamp time() const
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:69
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
double track1x
int hit_plane[kMaxHits]
void minuitFunctionz(int &nDim, double *gout, double &result, double par[], int flg)
double trkenddcosx[kMaxTrack]
void analyze(art::Event const &e) override
double translation
double z
int hit_channel[kMaxHits]
double Rho
Definition: doAna.cpp:13
double hit_ph[kMaxHits]
double track2y
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:89
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
double trkMCTruthT0[kMaxTrack]
RunNumber_t run() const
Definition: DataViewImpl.cc:82
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
Declaration of cluster object.
double trackgradientXZ[kMaxTrack]
double trketa_xy[kMaxTrack]
double trkendy[kMaxTrack]
Provides recob::Track data product.
double output[kMaxEvent][5][5]
Encapsulate the geometry of a wire.
double trkstarty[kMaxTrack]
const int kMaxTrack
p
Definition: test.py:223
double trkz[kMaxTrack][kMaxTrackHits]
art::ServiceHandle< geo::Geometry > fGeometry
Utility object to perform functions of association.
double trkPosx[kMaxTrack][3][1000]
std::string fMCTruthT0ModuleLabel
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
Encapsulate the construction of a single detector plane.
double trkd2[kMaxTrack]
int hit_trkid[kMaxHits]
double trkstartdcosy[kMaxTrack]
void minuitFunctionx(int &nDim, double *gout, double &result, double par[], int flg)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
tracking::Point_t Point_t
Definition: Track.h:53
double trkendz[kMaxTrack]
void endJob() override
double TrackPitchInView(recob::Track const &track, geo::View_t view, size_t trajectory_point=0U)
Returns the projected length of track on a wire pitch step [cm].
Definition: TrackUtils.cxx:78
object containing MC truth information necessary for making RawDigits and doing back tracking ...
double trkpitch[kMaxTrack][3]
double track1y
double trkstartz[kMaxTrack]
int trigger_offset(DetectorClocksData const &data)
static bool * b
Definition: config.cpp:1043
double trktheta[kMaxTrack]
double zdiff
double minx
double hit_peakT[kMaxHits]
EventNumber_t event() const
Definition: EventID.h:116
double zextrap(double y)
double trktheta_xz[kMaxTrack]
list x
Definition: train.py:276
std::string fTrkSpptAssocModuleLabel
double trkdQdxAverage[kMaxTrack]
double xdiff
int posminx
double trkdEdxSum[kMaxTrack]
TCEvent evt
Definition: DataStructs.cxx:7
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:291
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
double trkenddcosz[kMaxTrack]
double trklen[kMaxTrack]
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
Definition: TrackState.h:17
Utility functions to extract information from recob::Track
std::string fHitsModuleLabel
double trackgradientYZ[kMaxTrack]
double trkstartdcosz[kMaxTrack]
EventID id() const
Definition: Event.cc:37
double trkresrg[kMaxTrack][3][1000]
const int kMaxEvent
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
double trkx[kMaxTrack][kMaxTrackHits]
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
double ydiff
QTextStream & endl(QTextStream &s)
Event finding and building.
Encapsulate the construction of a single detector plane.
std::vector< double > possiblex
double trkstartdcosx[kMaxTrack]
double trkrange[kMaxTrack][3]
double trkplaneid[kMaxTrack][3][1000]
double trkpitchHit[kMaxTrack][3][1000]
int hit_wire[kMaxHits]