HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoTrackingPh.cxx
Go to the documentation of this file.
1 //
2 // File HarpoTrackingPh.cxx
3 //
12 #include "HarpoTrackingPh.h"
13 #include "HarpoConfig.h"
14 #include "HarpoDetSet.h"
15 #include "HarpoDebug.h"
16 #include "HarpoDccEvent.h"
17 #include "Pmm2Event.h"
18 #include "HarpoEvent.h"
19 #include "HarpoRecoEvent.h"
20 
21 #include "TFile.h"
22 #include "TStyle.h"
23 #include "TCanvas.h"
24 #include "TLatex.h"
25 #include "TGraph.h"
26 #include "TF1.h"
27 #include "TArrow.h"
28 #include "TLinearFitter.h"
29 #include "TMath.h"
30 #include "TSystem.h"
31 
32 #include <cstdlib>
33 #include <cstring>
34 #include <cassert>
35 #include <fstream>
36 #include <iostream>
37 
38 
39 ClassImp(HarpoTrackingPh);
40 
41 void HarpoTrackingPh::print()
42  {
43  unsigned int nd; // number of detectors
44  HarpoEventHeader *hdr = fEvt->GetHeader();
45 
46  assert(hdr != NULL);
47  hdr->print();
48 
49  for (nd = 0; nd < gkNDetectors; nd++) {
50  // if (fEvt->isdataExist(nd)) {
51  HarpoDccMap *plane = fEvt->GetDccMap(nd);
52  if (plane != NULL )plane->print();
53  }
54  }
55 
56 void HarpoTrackingPh::process()
57 {
58  // Bool_t drawEvent = kFALSE;
59  nEvents++;
60  if(gHarpoDebug>0)
61  Info("process","Event %ld",nEvents);
62 
63 
64 
65 
66 
67  for(Int_t i = 0; i< NTRACK; i++)
68  fNclTrack[i] = 0;
69  if(!fEvt->GetRecoEvent()) // Reconstructed data (clusters, tracks, matched tracks)
70  return;
71  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
72  reco->SetTrackType(KALMAN);
73  // Clusters in the event (both X and Y directions)
74  TClonesArray* clArray = reco->GetClustersArray();
75  // Int_t nCl = clArray->GetEntries();
76 
77 
78  // TMatrixD Xorig(4,1);
79  Double_t inc = fResKalman; // Space resolution
80  Double_t theta = fScatKalman; // Multiple scattering
81  TMatrixD Corig(4,4);
82  for(int a = 0;a<4;a++){
83  for(int b =0;b<4;b++){
84  if(a==b){
85  if(a<2)
86  Corig[a][b]=inc*inc; //initialisation de C
87  else
88  Corig[a][b]=(2*inc)*(2*inc)+theta*theta;
89  }
90  else
91  Corig[a][b]=0;
92  }
93  }
94 
95  for(Int_t plane = 0; plane<2; plane++){
96  if(!InitPlane(clArray,plane)) continue;
97 
98  Int_t counter = 0, color = 51, nClLeftOld = 0, nTry = 0;
99  while(1){
100  if(counter>50){
101  //Warning("process","Too many loops on plane %d",plane);
102  break;
103  }
104  if(gHarpoDebug>1)
105  Info("process","Look for seeds on plane %d %d",plane,counter);
106  counter++;
107 
108  Int_t iMin = NCHAN, iMax = 0, jMin = NADC, jMax = 0;
109  Int_t nClLeft = GetMapEdges(clArray,plane,iMin,iMax,jMin,jMax);
110  if(nClLeft<10){
111  if(gHarpoDebug>0)
112  Info("process","Only %d clusters left",nClLeft);
113  break;
114  }
115  if(nClLeft==nClLeftOld){
116  if(gHarpoDebug>0)
117  Info("process","No new tracks, %d clusters left",nClLeft);
118  break;
119  }
120  nClLeftOld = nClLeft;
121 
122 
123 
124 
125  for(Int_t side = 0; side<4; side++){ // Loop over starting map side
126  Int_t type = -1, index0 = 0, index1 = 0, imax0 = 0, imax1 = 0, sign = 0;
127  switch(side){
128  case 0: type = TCLUSTER; index0 = jMax; sign = -1; break;
129  case 1: type = CCLUSTER; index0 = iMax; sign = -1; break;
130  case 2: type = CCLUSTER; index0 = iMin; sign = +1; break;
131  case 3: type = TCLUSTER; index0 = jMin; sign = +1; break;
132  }
133  // if(side == 1 || side == 2) gHarpoDebug = 3;
134  // else continue;//gHarpoDebug = 0;
135  index1 = index0 + sign;
136  if(type == TCLUSTER){
137  if(index0<0 || index0>=NADC) continue;
138  while(NTcl[index1]<=0 && TMath::Abs(index1-index0) < 50) index1 += sign;
139  if(index1<0 || index1>=NADC) continue;
140  imax0 = NTcl[index0];
141  imax1 = NTcl[index1];
142  }else{
143  if(index0<0 || index0>=NCHAN) continue;
144  while(NCcl[index1]<=0 && TMath::Abs(index1-index0) < 50) index1 += sign;
145  if(index1<0 || index1>=NCHAN) continue;
146  imax0 = NCcl[index0];
147  imax1 = NCcl[index1];
148  }
149 
150 
151  if(TMath::Abs(index1-index0) >= 10){
152  for(Int_t i0 = 0; i0<imax0; i0++){
153  if(i0 == 20) break;
154  Int_t icl0;
155  if(type == CCLUSTER) icl0 = Ccl[index0][i0];
156  else icl0 = Tcl[index0][i0];
157  // used[icl0] = 3;
158  used[icl0] = -2000;
159  }
160  if(gHarpoDebug>1) Info("process","side %d: i0 = %d, i1 = %d",side,index0,index1);
161  continue;
162  }
163 
164 
165  if(gHarpoDebug>1)
166  Info("process","side %d => index = %d %d (%d %d)",side,index0,index1,imax0,imax1);
167  for(Int_t i0 = 0; i0<imax0; i0++){ // Loop over starting cluster pairs
168  if(i0 == 20) break;
169  Int_t icl0;
170  if(type == CCLUSTER) icl0 = Ccl[index0][i0];
171  else icl0 = Tcl[index0][i0];
172  if(gHarpoDebug>1) Info("process","used[%d] = %d",icl0,used[icl0]);
173  // if(used[icl0]) continue;
174  // used[icl0] = 3;
175  if(used[icl0] != -1000) continue;
176  used[icl0] = -2000;
177  for(Int_t i1 = 0; i1<imax1; i1++){
178  Int_t icl1;
179  if(type == CCLUSTER) icl1 = Ccl[index1][i1];
180  else icl1 = Tcl[index1][i1];
181  if(gHarpoDebug>1) Info("process","used[%d] = %d",icl1,used[icl1]);
182  // if(used[icl1]) continue;
183  if(used[icl1] != -1000) continue;
184  //used[icl1] = 3;
185 
186  //if(nTry==500){
187  // Warning("process","TOO MANY TRACK SEEDS");
188  // return;
189  //}
190  if(gHarpoDebug>1) Info("process","nTry = %d",nTry);
191  FindTrack(clArray,icl0,icl1,Corig,plane,color);
192  nTry++;
193 
194  }
195  }
196  }
197 
198  }
199  SpliceTracks(plane);
200  //SpliceTracksNew(plane);
201 
202  if(fHist[plane]){
203  fHist[plane]->GetXaxis()->SetRangeUser(0,511);
204  fHist[plane]->GetYaxis()->SetRangeUser(0,304);
205  fCanvas[plane]->Update();
206  }
207  }
208 }
209 
210 Int_t HarpoTrackingPh::InitPlane(TClonesArray* clArray, Int_t plane)
211 {
212 
213  fNtr[plane] = 0;
214 
215 
216  for(Int_t itr = 0; itr<NTRACK; itr++){
217  // if(fGraphs[plane][itr] != 0){
218  // fGraphs[plane][itr]->Delete();
219  // fGraphs[plane][itr] = new TGraph();
220  // }
221  fStartIndex[itr] = 0;
222  fEndIndex[itr] = 0;
223  fId[itr][0] = -1;
224  fId[itr][1] = -1;
225  for(Int_t j = 0; j<10;j++){
226  fStartPointX[itr][j] = -1000;
227  fStartPointZ[itr][j] = -1000;
228  // fStartAngle[itr][j] = -1000;
229  fStartDirX[itr][j] = -1000;
230  fStartDirZ[itr][j] = -1000;
231  // fEndPointX[itr][j] = -1000;
232  // fEndPointZ[itr][j] = -1000;
233  // fEndDirX[itr][j] = -1000;
234  // fEndDirZ[itr][j] = -1000;
235  // fEndAngle[itr][j] = -1000;
236  }
237  }
238  for(Int_t i = 0; i<NADC; i++){
239  NTcl[i] = 0;
240  for(Int_t k = 0; k<20; k++)
241  Tcl[i][k] = -1;
242  }
243  for(Int_t i = 0; i<NCHAN; i++){
244  NCcl[i] = 0;
245  for(Int_t k = 0; k<20; k++)
246  Ccl[i][k] = -1;
247  }
248 
249 
250  HarpoDccMap* map = fEvt->GetDccMap(plane);
251 
252  fQtot = 0;
253  fQused = 0;
254  for(Int_t i = 0; i<NCHAN; i++){
255  for(Int_t j = 0; j<NADC; j++){
256  // fMapTmp[i][j] = map->GetData(i,j);
257  fMapTmp[i][j] = 0;
258  if(map->GetData(i,j)>-500) fQtot += map->GetData(i,j);
259  }
260  }
261 
262  Int_t nCl = clArray->GetEntries();
263  if(nCl>=4000){
264  Warning("InitPlane","Too many clusters %d",nCl);
265  return 0;
266  }
267  // Organise Clusters in layers
268  for(Int_t icl = 0; icl<nCl; icl++){
269  reused[icl] = 0;
270  used[icl] = -1000;
271  fTrack[icl] = -1000;
272  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
273  Int_t index = cluster->GetIndex();
274  // Info("","index = %d",index);
275  //std::cout << plane << " " << cluster->GetPlane() << std::endl;
276  if(cluster->GetPlane() != plane) continue;
277  if(cluster->GetQ() < fQmin) continue;
278  if(cluster->GetType() == TCLUSTER && cluster->GetWidth() < fWidthMinT) continue;
279  if(cluster->GetType() == CCLUSTER && cluster->GetWidth() < fWidthMinC) continue;
280  if(cluster->GetWidth() > fWidthMax) continue;
281  if(cluster->GetType() == CCLUSTER){
282  if(NCcl[index]>=20) continue;
283  Ccl[index][NCcl[index]] = icl;
284  NCcl[index]++;
285  }
286  if(cluster->GetType() == TCLUSTER){
287  if(NTcl[index]>=20) continue;
288  Tcl[index][NTcl[index]] = icl;
289  NTcl[index]++;
290  }
291  }
292 
293  if(gHarpoDebug>2){
294  for(Int_t i = 0; i<NCHAN; i++){
295  if(NCcl[i]>0) Info("process","NCcl[%d] = %d",i,NCcl[i]);
296  }
297  for(Int_t i = 0; i<NADC; i++){
298  if(NTcl[i]>0) Info("process","NTcl[%d] = %d",i,NTcl[i]);
299  }
300  }
301 
302  return 1;
303 }
304 
305 
306 
307 Int_t HarpoTrackingPh::GetMapEdges(TClonesArray* clArray, Int_t plane, Int_t &iMin, Int_t &iMax, Int_t &jMin, Int_t &jMax)
308 {
309 
310  Int_t nClLeft = 0;
311  iMin = NCHAN+1; iMax = -1; jMin = NADC+1; jMax = -1;
312  Int_t nCl = clArray->GetEntries();
313  for(Int_t icl = 0; icl<nCl; icl++){
314  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
315  // if(cluster->GetIdClusterTrack(0)!=-1) continue;
316  // if(used[icl]) continue;
317  if(used[icl] != -1000) continue;
318  Int_t index = cluster->GetIndex();
319  if(cluster->GetQuality() != 0) continue;
320  if(cluster->GetPlane() != plane) continue;
321  if(cluster->GetQ() < fQmin) continue;
322  if(cluster->GetType() == TCLUSTER && cluster->GetWidth() < fWidthMinT) continue;
323  if(cluster->GetType() == CCLUSTER && cluster->GetWidth() < fWidthMinC) continue;
324  if(cluster->GetWidth() > fWidthMax) continue;
325  if(cluster->GetType() == TCLUSTER){
326  if(index>jMax) jMax = index;
327  if(index<jMin) jMin = index;
328  }
329  if(cluster->GetType() == CCLUSTER){
330  if(index>iMax) iMax = index;
331  if(index<iMin) iMin = index;
332  }
333  nClLeft++;
334  }
335 
336  if(gHarpoDebug>1)
337  Info("process","%d %d %d %d (%d)",iMin, iMax, jMin, jMax, nClLeft);
338 
339  return nClLeft;
340 }
341 
342 void HarpoTrackingPh::FindTrack(TClonesArray* clArray, Int_t icl0, Int_t icl1, TMatrixD Corig, Int_t plane, Int_t &color)
343 {
344 
345  if(gHarpoDebug>2)
346  Info("FindTrack","%d %d",icl0,icl1);
347  if(fNtr[plane]>=NTRACK) return;
348 
349 
350  if(fGraphs[plane][fNtr[plane]] != 0){
351  fGraphs[plane][fNtr[plane]]->Delete();
352  fGraphs[plane][fNtr[plane]] = new TGraph();
353  }
354 
355  HarpoRecoClusters* cluster0 = (HarpoRecoClusters*)clArray->At(icl0);
356  HarpoRecoClusters* cluster1 = (HarpoRecoClusters*)clArray->At(icl1);
357  Int_t type = cluster0->GetType();
358  if(type != 0 && type != 1) return;
359  // if(cluster0->GetIdClusterTrack(0)!=-1 && cluster1->GetIdClusterTrack(0)!=-1) continue;
360  Double_t x0 = cluster0->GetMean();
361  Double_t x1 = cluster1->GetMean();
362  Double_t index0 = cluster0->GetIndex();
363  Double_t index1 = cluster1->GetIndex();
364  Int_t q0 = cluster0->GetQ();
365  // Int_t itr = fNtr[plane];
366  if(TMath::Abs(x0-x1)>15) return;
367  TMatrixD Xorig(4,1);
368  Xorig[type][0] = x0;
369  Xorig[1-type][0] = index0;
370  Xorig[2+type][0] = x1-x0;
371  Xorig[3-type][0] = index1-index0;
372  if(gHarpoDebug>2)
373  Info("process","FindNext(%f,%f,%f,%f,%d)",Xorig[0][0],Xorig[1][0],Xorig[2][0],Xorig[3][0],fNtr[plane]);
374  if(TMath::Abs(Xorig[2+type][0]/Xorig[3-type][0])>1.5) return; // Angular cut for each cluster type
375 
376  if(fCanvas[plane]){
377  TArrow* arrow = new TArrow();
378  arrow->SetLineColor(kBlack);
379  arrow->SetFillColor(color);
380  fCanvas[plane]->cd();
381  arrow->DrawArrow(Xorig[0][0],Xorig[1][0],Xorig[0][0]+Xorig[2][0],Xorig[1][0]+Xorig[3][0],0.01);
382  fHist[plane]->GetXaxis()->SetRangeUser(Xorig[0][0] - 10,Xorig[0][0] + 10);
383  fHist[plane]->GetYaxis()->SetRangeUser(Xorig[1][0] - 10,Xorig[1][0] + 10);
384  fCanvas[plane]->Update();
385  }
386 
387 
388  // TArrayI* arr = new TArrayI(1000);
389  // arr->AddAt(icl0+1,0); // Keep 0 and 1 for start and end
390  // arr =
391  // FindNext(Xorig,Corig,fNtr[plane],clArray,arr, plane,color,3001,0);
392  // FindNext(Xorig,Corig,fNtr[plane],clArray,0, plane,color,3001,0);
393  Double_t angleorig = TMath::ATan2(Xorig[2][0],Xorig[3][0]);
394  FindNext(Xorig,Corig,angleorig,fNtr[plane],clArray,0,1,plane,color,3001,0,q0);
395  if(gHarpoDebug>1)
396  Info("process","Found new track %d",fNtr[plane]);
397  // Int_t nClTr = AddTrack(arr,clArray,plane);
398  Int_t nClTr = AddTrack(clArray,plane);
399  if(nClTr>10)
400  color += 4;
401  // delete arr;
402  if(gHarpoDebug>1)
403  Info("process","Track %d: %d clusters",fNtr[plane]-1,nClTr);
404  // fNtr[plane]++;
405 
406 
407 }
408 
409 
410 
411 
412 //Int_t HarpoTrackingPh::AddTrack(Int_t iTr, TClonesArray* clArray, Int_t plane)
413 //Int_t HarpoTrackingPh::AddTrack(TArrayI* array, TClonesArray* clArray, Int_t plane)
414 Int_t HarpoTrackingPh::AddTrack(TClonesArray* clArray, Int_t plane)
415 {
416 
417  Int_t n = 0;
418  if(fNtr[plane]<0 || fNtr[plane]>=NTRACK){
419  Warning("AddTrack","Wrong track number %d",fNtr[plane]);
420  return -10000;
421  }
422 
423  TMatrixD X(4,1);
424  Int_t index = (fStartIndex[fNtr[plane]]+9)%10;
425  //Info("AddTrack","%d %d %d",plane,fNtr[plane],index);
426  X[0][0] = fStartPointZ[fNtr[plane]][index] + fStartDirZ[fNtr[plane]][index];
427  X[1][0] = fStartPointX[fNtr[plane]][index] + fStartDirX[fNtr[plane]][index];
428  X[2][0] = -fStartDirZ[fNtr[plane]][index];
429  X[3][0] = -fStartDirX[fNtr[plane]][index];
430 
431  if(gHarpoDebug>1)
432  Info("AddTrack","New seed: (%g, %g, %g, %g)",X[0][0],X[1][0],X[2][0],X[3][0]);
433 
434  TMatrixD Corig(4,4);
435  for(int a = 0;a<4;a++){
436  for(int b =0;b<4;b++){
437  if(a==b){
438  if(a<2)
439  Corig[a][b]=fResKalman*fResKalman; //initialisation de C
440  else
441  Corig[a][b]=(2*fResKalman)*(2*fResKalman)+fScatKalman*fScatKalman;
442  }
443  else
444  Corig[a][b]=0;
445  }
446  }
447  TArrayI* arr = new TArrayI(2000);
448  Double_t angleorig = TMath::ATan2(X[2][0],X[3][0]);
449  arr = FindNextClosest(X, Corig, angleorig,fNtr[plane], clArray, arr, 1,plane, kGray, 3002, 2, -1);
450 
451  for(Int_t i = 0; i<2000; i++){
452  if(arr->At(i)-1<0) break;
453  n++;
454  }
455  //Info("AddTrack","Ncl = %d",n);
456  if(n<fNclMin){
457  delete arr;
458  return -n;
459  }
460 
461  if(fGraphs[plane][fNtr[plane]])
462  fGraphs[plane][fNtr[plane]]->Delete();
463  fGraphs[plane][fNtr[plane]] = new TGraph();
464  if(plane == XPLANE)
465  fGraphs[plane][fNtr[plane]]->SetName(Form("TrackX%d",fNtr[plane]));
466  else
467  fGraphs[plane][fNtr[plane]]->SetName(Form("TrackY%d",fNtr[plane]));
468 
469  for(Int_t i = 0; i<n; i++){
470  Int_t icl = arr->At(i) - 1;
471  //Info("AddTarack","icl = %d",icl);
472  if(icl<0) break;
473  if(icl>=4000) continue;
474  fTrack[icl] = fNtr[plane];
475  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
476  Double_t mean = cluster->GetMean();
477  Int_t index = cluster->GetIndex();
478  if(cluster->GetQuality() != 0) continue;
479  //cluster->AddIdClusterTrack(fNtr[plane]);
480  if(cluster->GetType() == TCLUSTER){
481  fGraphs[plane][fNtr[plane]]->SetPoint(fGraphs[plane][fNtr[plane]]->GetN(),index,mean);
482 
483  }
484  if(cluster->GetType() == CCLUSTER){
485  fGraphs[plane][fNtr[plane]]->SetPoint(fGraphs[plane][fNtr[plane]]->GetN(),mean,index);
486  }
487  }
488 
489  delete arr;
490 
491 
492  // if(gHarpoDebug>1)
493  // Info("AddTrack","Added track %d (plane %d) (%d cl) start(%g, %g) end(%g, %g)",fNtr[plane],plane,n,fStartPointX[fNtr[plane]][fStartIndex[fNtr[plane]]%10],fStartPointZ[fNtr[plane]][fStartIndex[fNtr[plane]]%10],fEndPointX[fNtr[plane]][fEndIndex[fNtr[plane]]%10],fEndPointZ[fNtr[plane]][fEndIndex[fNtr[plane]]%10]);
494 
495  fNtr[plane]++;
496 
497  return n;
498 }
499 
500 
501 
502 void HarpoTrackingPh::SpliceTracks(Int_t plane)
503 {
504  // return;
505 
506  //Double_t fDmin = 5, fThetaMin = -0.5;
507  Double_t fThetaMin = -0.5;
508  Double_t fDminForward = 20, fDminTransverse = 2;
509  // Double_t fDmin2 = fDmin*fDmin;
510  Int_t nSplice = 0;
511  Int_t fOffset = 0;
512  if(gHarpoDebug>1)
513  Info("SpliceTracks","plane %d: Ntracks = %d",plane,fNtr[plane]);
514 
515  Double_t x[fNtr[plane]][2], z[fNtr[plane]][2], th[fNtr[plane]][2];
516  Double_t px[fNtr[plane]][2], pz[fNtr[plane]][2];
517  Int_t spliced[fNtr[plane]];
518  GetQtracks(plane);
519  for(Int_t itr = 0; itr<fNtr[plane]; itr++){
520  spliced[itr] = itr;
521  if(fGraphs[plane][itr]->GetN()<=fOffset+1) continue;
522  fGraphs[plane][itr]->GetPoint(0 + fOffset,z[itr][0],x[itr][0]);
523  Double_t xTmp, zTmp;
524  fGraphs[plane][itr]->GetPoint(1 + fOffset,xTmp,zTmp);
525  px[itr][0] = xTmp - x[itr][0];
526  pz[itr][0] = zTmp - z[itr][0];
527  fGraphs[plane][itr]->GetPoint(fGraphs[plane][itr]->GetN()-1-fOffset,z[itr][1],x[itr][1]);
528  fGraphs[plane][itr]->GetPoint(fGraphs[plane][itr]->GetN()-2-fOffset,xTmp,zTmp);
529  px[itr][1] = xTmp - x[itr][1];
530  pz[itr][1] = zTmp - z[itr][1];
531  }
532  Int_t nLoops = 0;
533  while(1){
534  if(nLoops>10) break;
535  nLoops++;
536  Int_t nSplice = 0;
537 
538  for(Int_t itr1 = 0; itr1<fNtr[plane]; itr1++){
539  if(spliced[itr1] != itr1) continue;
540  for(Int_t itr2 = itr1+1; itr2<fNtr[plane]; itr2++){
541  if(spliced[itr2] != itr2) continue;
542  Double_t dMin = 1e8;
543  Int_t iMin = -1, jMin = -1;
544  for(Int_t i = 0; i<2; i++){
545  for(Int_t j = 0; j<2; j++){
546  Double_t d = TMath::Abs(x[itr1][i]-x[itr2][j]) + TMath::Abs(z[itr1][i]-z[itr2][j]);
547  hDistTracks->Fill(d);
548  if(d<dMin){
549  dMin = d;
550  iMin = i;
551  jMin = j;
552  }
553  }
554  }
555  Double_t dx1 = x[itr1][1-iMin] - x[itr1][iMin];
556  Double_t dz1 = z[itr1][1-iMin] - z[itr1][iMin];
557  Double_t dx2 = x[itr2][1-jMin] - x[itr2][jMin];
558  Double_t dz2 = z[itr2][1-jMin] - z[itr2][jMin];
559 
560 
561  Double_t cross = dx1*dx2 + dz1*dz2;
562  Double_t d1 = dx1*dx1 + dz1*dz1;
563  Double_t d2 = dx2*dx2 + dz2*dz2;
564  Double_t sd1 = TMath::Sqrt(d1);
565  Double_t sd2 = TMath::Sqrt(d2);
566  if(d1*d2 == 0) continue;
567  Double_t costheta = cross/TMath::Sqrt(d1*d2);
568 
569 
570  if(gHarpoDebug>1)
571  Info("SpliceTracks","%d %d: (d=%g,theta=%g)",itr1,itr2,dMin,costheta);
572  // if(iMin==jMin) costheta = -costheta;
573  if(costheta > fThetaMin) continue;
574 
575 
576  hDistMinTracks->Fill(dMin);
577  hThetaTracks->Fill(1-costheta);
578  hThetaTracksNtr->Fill(dMin,1-costheta);
579  Double_t dForward1 = TMath::Abs((x[itr1][iMin] - x[itr2][jMin])*dx1 +
580  (z[itr1][iMin] - z[itr2][jMin])*dz1)/sd1;
581  Double_t dTransverse1 = TMath::Abs((x[itr1][iMin] - x[itr2][jMin])*dz1 -
582  (z[itr1][iMin] - z[itr2][jMin])*dx1)/sd1;
583  Double_t dForward2 = TMath::Abs((x[itr1][iMin] - x[itr2][jMin])*dx2 +
584  (z[itr1][iMin] - z[itr2][jMin])*dz2)/sd2;
585  Double_t dTransverse2 = TMath::Abs((x[itr1][iMin] - x[itr2][jMin])*dz2 -
586  (z[itr1][iMin] - z[itr2][jMin])*dx2)/sd2;
587 
588  if(gHarpoDebug>1)
589  Info("SpliceTracks","%d %d: (dF=%g,dT=%g,theta=%g)",itr1,itr2,dForward1,dTransverse1,costheta);
590 
591  if(fCanvas[plane]){
592  Double_t xRect[4] = {x[itr1][iMin] + ( fDminForward*dx1/sd1 - fDminTransverse*dz1/sd1),
593  x[itr1][iMin] + (- fDminForward*dx1/sd1 - fDminTransverse*dz1/sd1),
594  x[itr1][iMin] + (- fDminForward*dx1/sd1 + fDminTransverse*dz1/sd1),
595  x[itr1][iMin] + ( fDminForward*dx1/sd1 + fDminTransverse*dz1/sd1)};
596  Double_t zRect[4] = {z[itr1][iMin] + ( fDminForward*dz1/sd1 + fDminTransverse*dx1/sd1),
597  z[itr1][iMin] + (- fDminForward*dz1/sd1 + fDminTransverse*dx1/sd1),
598  z[itr1][iMin] + (- fDminForward*dz1/sd1 - fDminTransverse*dx1/sd1),
599  z[itr1][iMin] + ( fDminForward*dz1/sd1 - fDminTransverse*dx1/sd1)};
600 
601  std::cout << "----------------" << std::endl;
602  for(Int_t i = 0; i<4; i++) std::cout << xRect[i] << " ";
603  std::cout << std::endl;
604  for(Int_t i = 0; i<4; i++) std::cout << zRect[i] << " ";
605  std::cout << std::endl;
606  TPolyLine* rect = new TPolyLine(4,zRect,xRect,"FL");
607  rect->SetFillColor(kRed);
608  rect->SetLineWidth(2);
609  fCanvas[plane]->cd();
610  rect->Draw("L");
611  fHist[plane]->GetYaxis()->SetRangeUser(x[itr1][iMin]-20,x[itr1][iMin]+20);
612  fHist[plane]->GetXaxis()->SetRangeUser(z[itr1][iMin]-20,z[itr1][iMin]+20);
613  fCanvas[plane]->Update();
614  gSystem->Sleep(100);
615  }
616 
617 
618  if(nLoops<10)
619  hThetaDist[nLoops]->Fill(dMin,1-costheta);
620  if(((dForward1<fDminForward && dTransverse1<fDminTransverse) ||
621  (dForward2<fDminForward && dTransverse2<fDminTransverse)) &&
622  iMin >= 0){
623 
624 
625  if(fCanvas[plane]){
626  //a1 = 1; a2 = 1;
627  fCanvas[plane]->cd();
628  TArrow* arrow1 = new TArrow();
629  arrow1->SetLineColor(kBlue);
630  arrow1->DrawArrow(z[itr1][iMin],x[itr1][iMin],z[itr1][iMin] + dz1,x[itr1][iMin]+dx1);
631  TArrow* arrow2 = new TArrow();
632  arrow2->SetLineColor(kRed);
633  arrow2->DrawArrow(x[itr2][jMin],x[itr2][jMin],z[itr2][jMin]+dz2,x[itr2][jMin]+dx2);
634  fCanvas[plane]->Update();
635  gSystem->Sleep(1000);
636  fHist[plane]->GetXaxis()->SetRangeUser(90,420);
637  fHist[plane]->GetYaxis()->SetRangeUser(0,288);
638  }
639 
640 
641 
642  spliced[itr2] = spliced[itr1];
643  x[itr1][iMin] = x[itr2][1-jMin];
644  z[itr1][iMin] = z[itr2][1-jMin];
645  px[itr1][iMin] = px[itr2][1-jMin];
646  pz[itr1][iMin] = pz[itr2][1-jMin];
647  th[itr1][iMin] = th[itr2][1-jMin];
648  nSplice++;
649  if(gHarpoDebug>1)
650  Info("SpliceTrack","Splicing %d %d (dF=%g,dT=%g,theta=%g)",itr1,itr2,dForward1,dTransverse1,costheta);
651  continue;
652  }
653 
654  Double_t a1 = ((x[itr1][iMin] - x[itr2][jMin])*dz2 - (z[itr1][iMin] - z[itr2][jMin])*dx2)/(dz1*dx2 - dz2*dx1);
655  Double_t a2 = ((x[itr1][iMin] - x[itr2][jMin])*dz1 - (z[itr1][iMin] - z[itr2][jMin])*dx1)/(dz1*dx2 - dz2*dx1);
656 
657  if(gHarpoDebug>1)
658  Info("SpliceTracks","%d %d: %g/%g %g/%g",itr1,itr2,a1,sd1,a2,sd2);
659 
660 
661 
662  if( TMath::Abs(a1 - 0.5) > 1) continue;
663  if( TMath::Abs(a2 - 0.5) > 1) continue;
664  if( (( a1*sd1 < -40 ) || (a1*sd1 > sd1 + 40)) &&
665  (( a2*sd2 < -40 ) || (a2*sd2 > sd2 + 40)))
666  continue;
667 
668  if(fCanvas[plane]){
669  //a1 = 1; a2 = 1;
670  fCanvas[plane]->cd();
671  TArrow* arrow1 = new TArrow();
672  arrow1->SetLineColor(kBlue+3);
673  arrow1->DrawArrow(x[itr1][iMin],z[itr1][iMin],x[itr1][iMin] + a1*dx1,z[itr1][iMin]+a1*dz1);
674  TArrow* arrow2 = new TArrow();
675  arrow2->SetLineColor(kRed+3);
676  arrow2->DrawArrow(x[itr2][jMin],z[itr2][jMin],x[itr2][jMin]+a2*dx2,z[itr2][jMin]+a2*dz2);
677  fCanvas[plane]->Update();
678  gSystem->Sleep(1000);
679  fHist[plane]->GetXaxis()->SetRangeUser(90,420);
680  fHist[plane]->GetYaxis()->SetRangeUser(0,288);
681  }
682 
683  Double_t dMax = a1*sd1;
684  Int_t i0 = iMin, tr0 = itr1;
685  Int_t i1 = iMin, tr1 = itr1;
686  if((1-a1)*sd1 > dMax){
687  dMax = (1-a1)*sd1;
688  i1 = i0;
689  tr1 = tr0;
690  i0 = 1-iMin;
691  tr0 = itr1;
692  }
693  if(a2*sd2 > dMax){
694  dMax = a2*sd2;
695  i1 = i0;
696  tr1 = tr0;
697  i0 = jMin;
698  tr0 = itr2;
699  }
700  if((1-a2)*sd2 > dMax){
701  dMax = (1-a2)*sd2;
702  i1 = i0;
703  tr1 = tr0;
704  i0 = 1-jMin;
705  tr0 = itr2;
706  }
707 
708  spliced[itr2] = spliced[itr1];
709  x[itr1][0] = x[tr0][i0];
710  z[itr1][0] = z[tr0][i0];
711  px[itr1][0] = px[tr0][i0];
712  pz[itr1][0] = pz[tr0][i0];
713  th[itr1][0] = th[tr0][i0];
714  x[itr1][1] = x[tr1][i1];
715  z[itr1][1] = z[tr1][i1];
716  px[itr1][1] = px[tr1][i1];
717  pz[itr1][1] = pz[tr1][i1];
718  th[itr1][1] = th[tr1][i1];
719  nSplice++;
720  if(gHarpoDebug>1)
721  Info("SpliceTrack","** Splicing %d %d",itr1,itr2);
722  continue;
723  }
724  }
725  if(gHarpoDebug>1)
726  Info("SpliceTracks","** %d spliced (%d)",nSplice,nLoops);
727  if(nSplice == 0) break;
728  }
729 
730 
731 
732 
733 
734 
735  Int_t newTrIndex[fNtr[plane]];
736  Int_t nClTr[fNtr[plane]];
737  Double_t qTr[fNtr[plane]];
738  Int_t nSpliced = 0;
739  for(Int_t itr = 0; itr<fNtr[plane]; itr++){
740  nClTr[itr] = 0;
741  qTr[itr] = 0;
742  if(spliced[itr] != itr){
743  Int_t itrNew = spliced[itr];
744  while(spliced[itrNew] != itrNew)
745  itrNew = spliced[itrNew];
746  spliced[itr] = itrNew;
747  Int_t n = fGraphs[plane][itr]->GetN();
748  for(Int_t k = 0; k<n; k++){
749  Double_t xTmp, zTmp;
750  fGraphs[plane][itr]->GetPoint(k,xTmp,zTmp);
751  fGraphs[plane][itrNew]->SetPoint(fGraphs[plane][itrNew]->GetN(),xTmp,zTmp);
752  // fGraphs[plane][itr]->SetPoint(k,0,0);
753  }
754  fGraphs[plane][itr]->Set(0);
755  // fGraphs[plane][itr]->Delete();
756  }else{
757  nSpliced++;
758  }
759  newTrIndex[itr] = -1;//nSpliced;
760  }
761 
762 
763 
764  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
765  // Clusters in the event (both X and Y directions)
766  TClonesArray* clArray = reco->GetClustersArray();
767  Int_t nCl = clArray->GetEntries();
768  if(gHarpoDebug>1)
769  Info("Splice","ntot = %d",nCl);
770  // Int_t nClTot = 0, nClTot2 = 0;
771  Int_t nFinal = 0;
772  for(Int_t icl = 0; icl<nCl; icl++){
773  Int_t oldTrack = fTrack[icl];
774  if(oldTrack<0 || oldTrack>=fNtr[plane]) continue;
775  Int_t splicedTrack = spliced[oldTrack];
776  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
777  cluster->AddIdClusterTrack(splicedTrack);
778  nClTr[splicedTrack]++;
779  // nClNew[newTrackId]++;
780  }
781  GetQtracks(plane);
782  for(Int_t itr1 = 0; itr1<fNtr[plane]; itr1++){
783  qTr[itr1] = fQcommon[itr1][itr1];
784  for(Int_t itr2 = itr1+1; itr2<fNtr[plane]; itr2++){
785  if(gHarpoDebug>1)
786  Info("SpliceTracks","%d => %d (%g %g %g)",itr2,spliced[itr1],fQcommon[itr1][itr1],fQcommon[itr2][itr2],fQcommon[itr1][itr2]);
787  if(!(fQcommon[itr1][itr2] > 0.4*fQcommon[itr2][itr2] || fQcommon[itr1][itr2] > 0.4*fQcommon[itr1][itr1] || fQcommon[itr1][itr2] > 20000))
788  continue;
789  if(gHarpoDebug>1)
790  Info("SpliceTracks","OK");
791  spliced[itr2] = spliced[itr1];
792  nClTr[spliced[itr1]] += nClTr[itr2];
793  qTr[spliced[itr1]] += fQcommon[itr2][itr2];
794  nSplice++;
795  }
796  }
797  for(Int_t itr1 = 0; itr1<fNtr[plane]; itr1++){
798  if(itr1 == spliced[itr1] && nClTr[itr1]>fNclMin2 && qTr[itr1]>20000){
799  newTrIndex[itr1] = nFinal;
800  nFinal++;
801  }else
802  newTrIndex[itr1] = -1;
803 
804  }
805 
806  for(Int_t icl = 0; icl<nCl; icl++){
807  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
808  if(cluster->GetPlane() != plane) continue;
809  // Int_t nTr = cluster->GetNtr();
810  cluster->RemoveAllClusterTrack();
811  Int_t oldTrack = fTrack[icl];
812  if(oldTrack<0 || oldTrack>=fNtr[plane]) continue;
813  Int_t sp = spliced[oldTrack];
814  if(sp<0) continue;
815  if(newTrIndex[sp]<0) continue;
816  cluster->AddIdClusterTrack(newTrIndex[sp]);
817  }
818 
819  for(Int_t itr = 0; itr<fNtr[plane]; itr++){
820  if(gHarpoDebug>1)
821  Info("Splice","nCl(%d) = %d",itr,nClTr[itr]);
822 
823  if(spliced[itr] != itr) continue;
824  if(newTrIndex[itr] == -1) continue;
825  if(gHarpoDebug>1)
826  Info("Splice","nCl(%d) = %d / %d",itr,nClTr[itr],fGraphs[plane][itr]->GetN());
827  Double_t x0 = -1000, slope = -1000, chi2 = -1000;
828  if(fGraphs[plane][itr]->GetN()>10){
829  TLinearFitter* linfit = new TLinearFitter(1,"pol1","D");
830  linfit->AssignData(fGraphs[plane][itr]->GetN(), 1, fGraphs[plane][itr]->GetX(), fGraphs[plane][itr]->GetY());
831  linfit->EvalRobust();
832  x0 = linfit->GetParameter(0);
833  slope = linfit->GetParameter(1);
834  chi2 = linfit->GetChisquare();
835  linfit->Delete();
836  }
837 
838  Double_t qTrTmp = GetQtrack(itr,plane);
839  fQused += qTrTmp;
840 
841  if(gHarpoDebug>1)
842  Info("SpliceTracks","Qleft = %g / %g",fQtot-fQused,fQtot);
843 
844  //if(qTr<50000) continue;
845 
846  //Int_t ttz = 0, typeTrack = 0, IdTrackMatching = 0;
847  Int_t typeTrack = 0, IdTrackMatching = 0;
848  Double_t xSt = x[itr][0], zSt = z[itr][0], thSt = th[itr][0], xEnd = x[itr][1], zEnd = z[itr][1], thEnd = th[itr][1];
849  Double_t pxSt = px[itr][0], pzSt = pz[itr][0], pxEnd = px[itr][1], pzEnd = pz[itr][1];
850  if(gHarpoDebug>1)
851  Info("SpliceTracks","new Track: %d (%d cl, q=%g) => (%g, %g, %g) (%g, %g, %g)", newTrIndex[itr],nClTr[itr], qTrTmp, xSt, zSt, thSt, xEnd, zEnd, thEnd);
852  HarpoRecoKalmanTracks* tr = new HarpoRecoKalmanTracks(newTrIndex[itr],Int_t(qTrTmp),slope,x0,chi2,typeTrack,IdTrackMatching,nClTr[itr],plane,xSt,zSt,pxSt,pzSt,xEnd,zEnd,pxEnd,pzEnd);
853  // Info("SpliceTracks","Add graph %p",fGraphs[plane][itr]);
854  fEvt->GetRecoEvent()->AddTracks(tr);
855  fId[newTrIndex[itr]][plane] = itr;
856  tr->Delete();
857  }
858  if(gHarpoDebug>1)
859  Info("SpliceTracks","Done");
860 
861 }
862 
863 
864 
865 Double_t HarpoTrackingPh::GetQtrack(Int_t itr, Int_t plane)
866 {
867  HarpoDccMap* map = fEvt->GetDccMap(plane);
868 
869  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
870  TClonesArray* clArray = reco->GetClustersArray();
871  Int_t nCl = clArray->GetEntries();
872  Double_t qTr = 0;
873  for(Int_t icl = 0; icl<nCl; icl++){
874  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
875  if(cluster->GetPlane() != plane) continue;
876  if(!cluster->CheckIdClusterTrack(itr)) continue;
877  Int_t type = cluster->GetType();
878  Int_t ind = cluster->GetIndex();
879  Int_t ind2 = cluster->GetIfirst();
880  Int_t width = cluster->GetWidth();
881  for(Int_t k = ind2; k< ind2 + width; k++){
882  Int_t i=-1 , j = -1; // undefined cell
883  if(type == CCLUSTER){ i = ind; j = k;}
884  if(type == TCLUSTER){ i = k; j = ind;}
885  if(i<0 || i>=NCHAN) continue;
886  if(j<0 || j>=NADC) continue;
887  Double_t q = map->GetData(i,j);
888  // Double_t q = fMapTmp[i][j];
889  if(q<-500) continue;
890  ULong64_t test = 1 << itr;
891  //Info("GetQtrack","TEST = %llX | %llX => %llX, %llX",test,fMapTmp[i][j],fMapTmp[i][j]&test,fMapTmp[i][j]|test);
892  if(fMapTmp[i][j] & test) continue;
893  qTr += q;
894  // fMapTmp[i][j] = -1000;
895  fMapTmp[i][j] |= test;
896  }
897  }
898  if(gHarpoDebug>0)
899  Info("GetQtrack","Plane %d, Track %d: Q = %g/%g",plane,itr,qTr,fQtot);
900  return qTr;
901 }
902 
903 Double_t HarpoTrackingPh::GetQtracks(Int_t plane)
904 {
905  HarpoDccMap* map = fEvt->GetDccMap(plane);
906 
907  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
908  TClonesArray* clArray = reco->GetClustersArray();
909  Int_t nCl = clArray->GetEntries();
910  // Double_t qTr = 0;
911 
912  ULong_t used4Q[NCHAN][NADC];
913  for(Int_t i = 0; i<NCHAN; i++){
914  for(Int_t j = 0; j<NADC; j++){
915  // fMapTmp[i][j] = map->GetData(i,j);
916  used4Q[i][j] = 0;
917  }
918  }
919 
920  for(Int_t icl = 0; icl<nCl; icl++){
921  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
922  if(cluster->GetPlane() != plane) continue;
923  Int_t type = cluster->GetType();
924  Int_t ind = cluster->GetIndex();
925  Int_t ind2 = cluster->GetIfirst();
926  Int_t width = cluster->GetWidth();
927  Int_t nTr = cluster->GetNtr();
928  for(Int_t iTr = 0; iTr<nTr;iTr++){
929  Int_t itr = cluster->GetIdClusterTrack(iTr);
930  for(Int_t k = ind2; k< ind2 + width; k++){
931  Int_t i=-1 , j = -1; // undefined cell
932  if(type == CCLUSTER){ i = ind; j = k;}
933  if(type == TCLUSTER){ i = k; j = ind;}
934  if(i<0 || i>=NCHAN) continue;
935  if(j<0 || j>=NADC) continue;
936  Double_t q = map->GetData(i,j);
937  // Double_t q = fMapTmp[i][j];
938  if(q<-500) continue;
939  ULong64_t test = 1 << itr;
940  //Info("GetQtrack","TEST = %llX | %llX => %llX, %llX",test,fMapTmp[i][j],fMapTmp[i][j]&test,fMapTmp[i][j]|test);
941  if(used4Q[i][j] & test) continue;
942  // qTr += q;
943  // fMapTmp[i][j] = -1000;
944  used4Q[i][j] |= test;
945  }
946  }
947  }
948 
949  // Double_t qTr[fNtr[plane]][fNtr[plane]];
950  for(Int_t i = 0; i<fNtr[plane]; i++){
951  for(Int_t j = 0; j<fNtr[plane]; j++){
952  fQcommon[i][j] = 0;
953  }
954  }
955  for(Int_t i = 0; i<NCHAN; i++){
956  for(Int_t j = 0; j<NADC; j++){
957  Double_t q = map->GetData(i,j);
958  if(q<-500) continue;
959 
960  for(Int_t itr1 = 0; itr1<fNtr[plane]; itr1++){
961  ULong64_t test1 = 1 << itr1;
962  if(!(used4Q[i][j] & test1)) continue;
963  for(Int_t itr2 = itr1; itr2<fNtr[plane]; itr2++){
964  ULong64_t test2 = 1 << itr2;
965  if(!(used4Q[i][j] & test2)) continue;
966  fQcommon[itr1][itr2] += q;
967  }
968  }
969  }
970  }
971 
972  if(gHarpoDebug>1){
973  for(Int_t i = 0; i<fNtr[plane]; i++){
974  for(Int_t j = 0; j<fNtr[plane]; j++){
975  std::cout << fQcommon[i][j] << " ";
976  }
977  std::cout << std::endl;
978  }
979  }
980 
981  // Info("GetQtrack","Plane %d, Track %d: Q = %g/%g",plane,itr,qTr,fQtot);
982  return 0;
983 }
984 
985 TArrayI* HarpoTrackingPh::FindNext(TMatrixD X, TMatrixD C, Double_t angle, Int_t iTr, TClonesArray* clArray, TArrayI* arr, Int_t ncl, Int_t plane, Int_t color, Int_t fill, Int_t smooth, Int_t qold)
986 {
987  // return FindNextClosest(X,C,iTr,clArray,arr,plane,color,fill,smooth);
988  return FindNextClosest(X,C,angle,iTr,clArray,arr,ncl,plane,color,fill,smooth,qold);
989 }
990 
991 
992 TArrayI* HarpoTrackingPh::FindNextClosest(TMatrixD X, TMatrixD C, Double_t angle, Int_t iTr, TClonesArray* clArray, TArrayI* arr, Int_t ncl, Int_t plane, Int_t color, Int_t fill, Int_t smooth, Int_t qold)
993 {
994  if(gHarpoDebug>2) Info("FindNextClosest","Test2");
995  if(fCanvas[plane]) gSystem->Sleep(50);
996  if(gHarpoDebug>2) Info("FindNextClosest","angle = %g / %d = %g",angle,ncl,angle/ncl);
997 
998  TMatrixD F(4,4);
999  TMatrixD Xproj(4,1);
1000  // Double_t fNclMin = 20;
1001  Double_t deltaRes = fResFinder;
1002  Double_t deltaScat = fScatFinder;
1003  Double_t theta = fScatKalman; // Multiple scattering
1004  Int_t fNskipRows = 6;//10;
1005  Double_t fRangeMax = 5;
1006 
1007  if(gHarpoDebug>2)
1008  Info("FindNextClosest","%f %f %f %f",X[0][0],X[1][0],X[2][0],X[3][0]);
1009 
1010  Double_t angleOld = TMath::ATan2(X[2][0], X[3][0]);
1011  Int_t ntr = 0;
1012  Double_t distMin = 1000, meanMin = -10000, rangeMin = 0, angleMin = -1000, resMin = -1;
1013  Int_t iclMin = -1, typeMin = -1, iMin = -1, kMin = -1, qMin = -1;
1014  //Info("FindNextClosest","###############################");
1015  for(Int_t k = 1; k<=fNskipRows; k++){
1016  if(k>distMin + kMin) break;
1017  for(Int_t type = 0; type<2; type++){ // look for T- and C-clusters
1018  //if(gHarpoDebug>2)
1019  // if(TMath::Abs(X[2+type][0]/X[3-type][0])>fMaxSlopeType) continue; // Angular cut for each cluster type
1020  if(smooth != 2 && TMath::Abs(X[2+type][0]/X[3-type][0])>fMaxSlopeType) continue; // Angular cut for each cluster type
1021  // if(gHarpoDebug>2)
1022  // Info("FindNextClosest","type: %d",type);
1023 
1024  // Define next layer
1025  Int_t I;
1026  if(X[3-type][0]>0) I = Int_t(X[1-type][0]+0.5) + k;
1027  else I = Int_t(X[1-type][0]+0.5) - k;
1028  //Info("FindNextClosest","type: %d, I = %d",type,I);
1029  if(I<0 || (type==CCLUSTER && I>=NCHAN) || (type==TCLUSTER && I>=NADC)) continue;
1030  Int_t n = 0;
1031  if(type == TCLUSTER) n = NTcl[I];
1032  else n = NCcl[I];
1033  // if(gHarpoDebug>2)
1034  // Info("FindNextClosest","%d %d",I,n);
1035 
1036  for(int i = 0;i<4;i++){
1037  for(int j =0;j<4;j++){
1038  F[i][j]=0;
1039  if(i==j)
1040  F[i][j]=1; //initialisation de F
1041  }
1042  }
1043  F[0][2]=(I-X[1-type][0])*1./X[3-type][0];
1044  F[1][3]=(I-X[1-type][0])*1./X[3-type][0];
1045  Xproj=F*X;
1046 
1047  for(Int_t i = 0; i<n; i++){ // For all clusters in the layer
1048  if(i == 20) break;
1049  Int_t icl;
1050  if(type == TCLUSTER) icl = Tcl[I][i];
1051  else icl = Ccl[I][i];
1052  if(smooth == 0 && used[icl] > -1000) continue; // first pass: If already in a track, skip
1053  if(smooth == 1 && used[icl] > -1000 && reused[icl] != 0) continue; // second pass: If already in a track (except this one), skip
1054  if(smooth == 2 && used[icl] != iTr) continue; // reordering: If not in this track, skip
1055 
1056  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
1057  if(cluster->GetPlane() != plane) continue;
1058  if((cluster->GetQuality() & 3) == 3) continue;
1059  if(cluster->GetType() != type) Warning("FindNextClosest","??? %d != %d",cluster->GetType(), type);
1060  Double_t mean = cluster->GetMean();
1061  // Int_t index = cluster->GetIndex();
1062  Int_t q = cluster->GetQ();
1063 
1064  if(gHarpoDebug>2)
1065  Info("FindNextClosest","#proj: %f %f %f %f",Xproj[0][0],Xproj[1][0],Xproj[2][0],Xproj[3][0]);
1066 
1067 
1068  Double_t res = GetResolution(cluster);
1069  Double_t range = deltaRes*res + deltaScat*theta*k;
1070  if(range>fRangeMax) range = fRangeMax;
1071 
1072  Double_t dist = TMath::Abs(mean - Xproj[type][0]);// + k;
1073  //fHistQQtestVsDist->Fill(q,qold);
1074  fHistQQtestVsDist->Fill(q*1./qold,dist);
1075  Double_t angleTmp;
1076  if(type == 0) angleTmp = TMath::ATan2(mean - X[0][0], I - X[1][0] );
1077  else angleTmp = TMath::ATan2(I - X[0][0], mean - X[1][0]);
1078  Double_t dA = fmod(angleTmp - angle/ncl + 4*TMath::Pi(),2*TMath::Pi()) - TMath::Pi();
1079  if(TMath::Abs(dA)<TMath::Pi()/2 && !(smooth == 1 && used[icl] == iTr)) continue;
1080 
1081  Double_t dA2 = fmod(angleTmp - angleOld + 4*TMath::Pi(),2*TMath::Pi()) - TMath::Pi();
1082  //Info("","dA = %g - %g = %g",angleTmp,angle/ncl,dA);
1083  if(TMath::Abs(dA2)<TMath::Pi()*0.7 && !(smooth == 1 && used[icl] == iTr)) continue;
1084  if(smooth == 1 && used[icl] == iTr) dist = 0; // Catch
1085  //Info("FindNextClosest","** %d %g %g %g %d %d",icl, dist,mean,range,type,I);
1086  if((dist < range || smooth == 2) && (dist+k<distMin || (dist == 0 && k<=kMin))){
1087  //if((dist < range) && (dist+k<distMin || (dist == 0 && k<=kMin))){
1088  distMin = dist;
1089  iclMin = icl;
1090  meanMin = mean;
1091  rangeMin = range;
1092  typeMin = type;
1093  iMin = I;
1094  kMin = k;
1095  if(smooth == 1) kMin = 100;
1096  angleMin = angleTmp;
1097  resMin = res;
1098  qMin = q;
1099  }
1100  }
1101  }
1102  }
1103 
1104  fHistQQmin->Fill(qMin*1./qold,distMin);
1105 
1106  // Info("FindNextClosest","%d: %g %d %g %g %d %d %d %g %g",smooth,distMin,iclMin,meanMin,rangeMin,typeMin,iMin,kMin,angleMin,resMin);
1107  // Info("FindNextClosest","used[%d] = %d",iclMin,used[iclMin]);
1108 
1109  if(iclMin!=-1) ntr++;
1110  else{ // END OF TRACK
1111  if(smooth == 0 && fStartIndex[iTr]>10){ // end of first pass
1112  // X[2][0] = -X[2][0];
1113  // X[3][0] = -X[3][0];
1114 
1115  Int_t index = (fStartIndex[iTr] + 4)%10;
1116  X[0][0] = fStartPointZ[iTr][index];
1117  X[1][0] = fStartPointX[iTr][index];
1118  X[2][0] = -fStartDirZ[iTr][index];
1119  X[3][0] = -fStartDirX[iTr][index];
1120 
1121  if(gHarpoDebug>1)
1122  Info("FindNextClosest","End of first pass (Ncl = %d). New seed: (%g, %g, %g, %g)",ncl,X[0][0],X[1][0],X[2][0],X[3][0]);
1123 
1124  TMatrixD Corig(4,4);
1125  for(int a = 0;a<4;a++){
1126  for(int b =0;b<4;b++){
1127  if(a==b){
1128  if(a<2)
1129  Corig[a][b]=fResKalman*fResKalman; //initialisation de C
1130  else
1131  Corig[a][b]=(2*fResKalman)*(2*fResKalman)+fScatKalman*fScatKalman;
1132  }
1133  else
1134  Corig[a][b]=0;
1135  }
1136  }
1137  // ncl++;
1138  // Double_t angleorig = TMath::ATan2(X[2][0],X[3][0]);
1139  // arr = FindNextClosest(X, Corig, angleorig, iTr, clArray, arr, 1, plane, color, fill, 1);
1140  arr = FindNextClosest(X, Corig, -angle, iTr, clArray, arr, ncl, plane, color, fill, 1, qMin);
1141  }else{ // end of reverse pass
1142  // Int_t Ncl = 0;
1143  if(gHarpoDebug>1)
1144  Info("FindNextClosest","End of last pass (%d clusters)",ncl);
1145  // while(arr->At(Ncl) != 0){
1146  // Ncl++;
1147  // }
1148  // std::cout << "*** " << iTr << " => " << Ncl << " ***" << std::endl;
1149  }
1150  return arr;
1151  }
1152 
1153  if(gHarpoDebug>2)
1154  Info("FindNextClosest","Add New Cluster");
1155 
1156 
1157  TMatrixD H(2,4);
1158  TMatrixD Ht(2,4);
1159  TMatrixD Q(4,4);
1160  TMatrixD G(2,2);
1161  TMatrixD V(2,2);
1162  // TMatrixD Id(4,4);
1163 
1164 
1165  Double_t inc = resMin;//fResKalman; // Space resolution
1166  if(inc==0) inc = 0.000001;
1167  for(int i = 0;i<4;i++){
1168  for(int j =0;j<4;j++){
1169  if(i==j){
1170  if(i>1)
1171  Q[i][j]=theta*theta; //initialisation de Q
1172  else
1173  Q[i][j]=0;
1174  }
1175  }
1176  }
1177  for(int i = 0;i<2;i++){
1178  for(int j =0;j<2;j++){
1179  if(i==j){
1180  H[i][j]=1;
1181  V[i][j]=inc*inc;
1182  G[i][j]=1./inc*inc;
1183  }else{
1184  H[i][j]=0;
1185  V[i][j]=0;
1186  G[i][j]=0;
1187  }
1188  }
1189  }
1190  Ht=H;
1191  Ht.T();
1192  // G=V;
1193  // G.Invert();
1194 
1195 
1196  TMatrixD Ft(4,4);
1197  TMatrixD Cinv(4,4);
1198  TMatrixD M(2,1);
1199  TMatrixD Cproj(4,4);
1200  TMatrixD Xnew(4,1);
1201  TMatrixD Cnew(4,4);
1202  TMatrixD Cprojinv(4,4);
1203 
1204  for(int i = 0;i<4;i++){
1205  for(int j =0;j<4;j++){
1206  F[i][j]=0;
1207  if(i==j)
1208  F[i][j]=1; //initialisation de F
1209  }
1210  }
1211  F[0][2]=(iMin-X[1-typeMin][0])*1./X[3-typeMin][0];
1212  F[1][3]=(iMin-X[1-typeMin][0])*1./X[3-typeMin][0];
1213  Xproj=F*X;
1214  // Kalman matrices calculations
1215  Ft = F;
1216  Ft.T();
1217  M[1-typeMin][0] = iMin;
1218  M[typeMin][0] = meanMin;
1219  Cproj=F*C*Ft+Q;
1220  Cprojinv=Cproj;
1221  Cprojinv.Invert();
1222  Cinv=Cprojinv+Ht*G*H;
1223  Cnew=Cinv;
1224  Cnew.Invert();
1225  Xnew=Cnew*((Cprojinv*Xproj)+(Ht*G*M));
1226 
1227 
1228  if(fCanvas[plane]){ // Display for Reco monitor
1229  Float_t x[3] = {(Float_t) (X[1-typeMin][0]),
1230  (Float_t) (iMin),
1231  (Float_t) (iMin)};
1232  Float_t y[3] = {(Float_t) (X[typeMin][0]),
1233  (Float_t) (Xproj[typeMin][0]-rangeMin),
1234  (Float_t) (Xproj[typeMin][0]+rangeMin)};
1235  if(gHarpoDebug>2)
1236  Info("FindNextClosest","%g %g %g %g %g %g",x[0],y[0],x[1],y[1],x[2],y[2]);
1237  TPolyLine* triangle;
1238  if(typeMin == TCLUSTER) triangle = new TPolyLine(3,x,y,"FL");
1239  else triangle = new TPolyLine(3,y,x,"FL");
1240  triangle->SetLineColor(kBlack);
1241  //triangle->SetFillColor(kRed);
1242  triangle->SetFillColor(color);
1243  triangle->SetFillStyle(fill);
1244  fCanvas[plane]->cd();
1245  triangle->Draw("FL");
1246  fHist[plane]->GetXaxis()->SetRangeUser(X[0][0]-20,X[0][0]+20);
1247  fHist[plane]->GetYaxis()->SetRangeUser(X[1][0]-20,X[1][0]+20);
1248  fCanvas[plane]->Update();
1249  }
1250  used[iclMin] = iTr;
1251  if(smooth == 1) reused[iclMin] = 1;
1252  else reused[iclMin] = 0;
1253  if(smooth == 2) used[iclMin] = 1000;
1254 
1255  Bool_t skip = kFALSE; // If already in this track, skip
1256  Int_t kcl = 0;
1257  if(arr != 0){
1258  while(arr->At(kcl) != 0){
1259  if(arr->At(kcl) == 0) break;
1260  if(arr->At(kcl) == iclMin+1){
1261  skip = kTRUE;
1262  break;
1263  }
1264  kcl++;
1265  }
1266  if(!skip)
1267  //Info("FindNext","Add cluster %d",iclMin);
1268  arr->AddAt(iclMin+1,kcl);
1269  }
1270  if(gHarpoDebug>2)
1271  Info("FindNextClosest","Add cluster %d, (%d %p) %d",iclMin,iTr,arr,fStartIndex[iTr]%10);
1272 
1273  if(smooth == 0){
1274  fStartPointZ[iTr][fStartIndex[iTr]%10] = X[0][0];
1275  fStartPointX[iTr][fStartIndex[iTr]%10] = X[1][0];
1276  fStartDirZ[iTr][fStartIndex[iTr]%10] = X[2][0];
1277  fStartDirX[iTr][fStartIndex[iTr]%10] = X[3][0];
1278  fStartIndex[iTr]++;
1279  }//else{
1280  // fEndPointZ[iTr][fEndIndex[iTr]%10] = X[0][0];
1281  // fEndPointX[iTr][fEndIndex[iTr]%10] = X[1][0];
1282  // fEndDirZ[iTr][fStartIndex[iTr]%10] = X[2][0];
1283  // fEndDirX[iTr][fStartIndex[iTr]%10] = X[3][0];
1284  // fEndIndex[iTr]++;
1285  //}
1286 
1287  if(smooth == 2){
1288  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(iclMin);
1289  cluster->SetZfit(0,Xnew[0][0]);
1290  cluster->SetXfit(0,Xnew[1][0]);
1291  }
1292  angle += angleMin;
1293  ncl++;
1294  if(gHarpoDebug>2)
1295  Info("FindNextClosest","Test %g %d %p %p %d %d %d %d %d", angle, iTr, clArray, arr, ncl, plane, color, fill,smooth);
1296  arr = FindNextClosest(Xnew, Cnew, angle, iTr, clArray, arr, ncl, plane, color, fill,smooth,qMin);
1297 
1298 
1299  return arr;
1300 }
1301 
1302 
1303 Double_t HarpoTrackingPh::GetResolution(HarpoRecoClusters* cl)
1304 {
1305 
1306  if((cl->GetQuality() & 3) == 3) return 2;
1307  if(cl->GetQuality() & 4) return 1;
1308  if(cl->GetQuality() & 3) return 1;//0.5;
1309 
1310 
1311  if(cl->GetType() == CCLUSTER) return 1;
1312  if(cl->GetType() == TCLUSTER) return 0.2;
1313 
1314  return 10;
1315 }
1316 
1317 
1318  void HarpoTrackingPh::Init()
1319  {
1320 
1321  // Initialise histograms here
1322 
1323  fCanvas[0] = NULL;
1324  fCanvas[1] = NULL;
1325  fHist[0] = NULL;
1326  fHist[1] = NULL;
1327 
1328  fNclMin = 10;
1329  fNclMin2 = 15;
1330  fMaxSlopeType = 1.5;
1331  fResKalman = 0.2;
1332  fScatKalman = 1;
1333  fResFinder = 1;
1334  fScatFinder = 0.5;
1335 
1336  fQmin = 200;
1337  fWidthMinT = 2;
1338  fWidthMinC = 5;
1339  fWidthMax = 30;
1340 
1341 
1342 
1343 
1344  // fNtr = 0;
1345  hNcl = new TH1F("hNcl","hNcl",500,0,500);
1346 
1347  for(Int_t itr = 0; itr<NTRACK; itr++){
1348  fGraphs[0][itr] = 0;
1349  fGraphs[1][itr] = 0;
1350  }
1351 
1352 
1353  const Int_t nbinsDist = 200;
1354  Double_t xminDist = 1;
1355  Double_t xmaxDist = 10000;
1356  Double_t logxminDist = TMath::Log(xminDist);
1357  Double_t logxmaxDist = TMath::Log(xmaxDist);
1358  Double_t binwidthDist = (logxmaxDist-logxminDist)/nbinsDist;
1359  Double_t xbinsDist[nbinsDist+1];
1360  xbinsDist[0] = xminDist;
1361  for (Int_t i=1;i<=nbinsDist;i++)
1362  xbinsDist[i] = TMath::Exp(logxminDist+i*binwidthDist);
1363 
1364 
1365  const Int_t nbinsTheta = 100;
1366  Double_t xminTheta = 1e-5;
1367  Double_t xmaxTheta = 2;
1368  Double_t logxminTheta = TMath::Log(xminTheta);
1369  Double_t logxmaxTheta = TMath::Log(xmaxTheta);
1370  Double_t binwidthTheta = (logxmaxTheta-logxminTheta)/nbinsTheta;
1371  Double_t xbinsTheta[nbinsTheta+1];
1372  xbinsTheta[0] = xminTheta;
1373  for (Int_t i=1;i<=nbinsTheta;i++)
1374  xbinsTheta[i] = TMath::Exp(logxminTheta+i*binwidthTheta);
1375 
1376  hDistTracks = new TH1F("hDistTracks","",nbinsDist,xbinsDist);
1377  hDistMinTracks = new TH1F("hDistMinTracks","",nbinsDist,xbinsDist);
1378  hThetaTracks = new TH1F("hThetaTracks","",nbinsTheta,xbinsTheta);//,100,-1,1);
1379  // hThetaTracksNtr = new TH2F("hThetaTracksNtr","",NTRACK,0,NTRACK,100,-1,1);
1380  hThetaTracksNtr = new TH2F("hThetaTracksNtr","",nbinsDist,xbinsDist,nbinsTheta,xbinsTheta);
1381  for(Int_t i = 0; i<10; i++)
1382  hThetaDist[i] = new TH2F(Form("hThetaDist%d",i),"",nbinsDist,xbinsDist,nbinsTheta,xbinsTheta);
1383 
1384 
1385  const Int_t nbinsQ = 100;
1386  Double_t xminQ = 1e-3;
1387  Double_t xmaxQ = 1e3;
1388  Double_t logxminQ = TMath::Log(xminQ);
1389  Double_t logxmaxQ = TMath::Log(xmaxQ);
1390  Double_t binwidthQ = (logxmaxQ-logxminQ)/nbinsQ;
1391  Double_t xbinsQ[nbinsQ+1];
1392  xbinsQ[0] = xminQ;
1393  for (Int_t i=1;i<=nbinsQ;i++)
1394  xbinsQ[i] = TMath::Exp(logxminQ+i*binwidthQ);
1395  fHistQQtestVsDist = new TH2F("fHistQQtestVsDist","",nbinsQ,xbinsQ,nbinsDist,xbinsDist);
1396  fHistQQmin = new TH2F("fHistQQmin","",nbinsQ,xbinsQ,nbinsDist,xbinsDist);
1397 
1398  }
1399 
1400 void HarpoTrackingPh::Save(char * /* mode */)
1401  {
1402 
1403  // TString * hstFile = gHConfig->GetHistFile();
1404  // if ( hstFile == NULL ) {
1405  // std::cout << gHConfig->GetProgramName() << " " <<
1406  // "No Hist File name given, use default" << std::endl;
1407  // // const char* dir = gSystem->ExpandPathName("$HARPO_ANA_DIR/reco");
1408  // const char* dir = "$HARPO_ANA_DIR/reco";
1409  // if(gSystem->AccessPathName(dir))
1410  // hstFile = new TString(Form("trackingKalman%lli.root",gHConfig->GetRunNo()));
1411  // else
1412  // hstFile = new TString(Form("%s/trackingKalman%lli.root",dir,gHConfig->GetRunNo()));
1413  // // hstFile = new TString(Form("outputs/anaph%lli.root",gHConfig->GetRunNo()));
1414  // }
1415 
1416  // TFile *hf = new TFile(hstFile->Data(),"RECREATE");
1417  OpenHistFile("trackingKalman");
1418  // Save histograms here
1419  hDistTracks->Write();
1420  hDistMinTracks->Write();
1421  hThetaTracks->Write();
1422  hThetaTracksNtr->Write();
1423  for(Int_t i = 0; i<10; i++){
1424  if(hThetaDist[i])
1425  hThetaDist[i]->Write();
1426  }
1427  fHistQQtestVsDist->Write();
1428  fHistQQmin->Write();
1429 
1430  // hf->Close();
1431  // printf("fNewFile %s closed\n", hstFile->Data() );
1432 
1433  }
1434 
NTRACK
#define NTRACK
Definition: HarpoHoughTracking.h:22
HarpoTrackingPh::FindTrack
void FindTrack(TClonesArray *clArray, Int_t icl0, Int_t icl1, TMatrixD Corig, Int_t plane, Int_t &color)
Definition: HarpoTrackingPh.cxx:342
HarpoTrackingPh::fCanvas
TVirtualPad * fCanvas[2]
Definition: HarpoTrackingPh.h:110
HarpoEvent::GetRecoEvent
HarpoRecoEvent * GetRecoEvent()
Definition: HarpoEvent.h:46
HarpoTrackingPh::Save
void Save(char *mode=NULL)
Definition: HarpoTrackingPh.cxx:1400
XPLANE
#define XPLANE
Definition: HarpoConfig.h:24
NCHAN
#define NCHAN
Definition: HarpoDccMap.h:16
HarpoTrackingPh::fStartDirX
Double_t fStartDirX[NTRACK][10]
Definition: HarpoTrackingPh.h:90
HarpoTrackingPh::fId
Int_t fId[NTRACK][2]
Definition: HarpoTrackingPh.h:99
HarpoRecoClusters::GetMean
Double_t GetMean()
Definition: HarpoRecoEvent.h:39
HarpoTrackingPh::SpliceTracks
void SpliceTracks(Int_t plane)
Definition: HarpoTrackingPh.cxx:502
HarpoRecoClusters::AddIdClusterTrack
void AddIdClusterTrack(Int_t val)
Definition: HarpoRecoEvent.cxx:75
HarpoTrackingPh::NCcl
Int_t NCcl[NCHAN]
Definition: HarpoTrackingPh.h:78
HarpoDccMap::GetData
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
Definition: HarpoDccMap.cxx:84
HarpoTrackingPh::fStartPointZ
Double_t fStartPointZ[NTRACK][10]
Definition: HarpoTrackingPh.h:88
HarpoRecoEvent
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
Definition: HarpoRecoEvent.h:474
HarpoRecoEvent::AddTracks
void AddTracks(HarpoRecoTracks *val)
Definition: HarpoRecoEvent.h:526
HarpoTrackingPh::fStartIndex
Int_t fStartIndex[NTRACK]
Definition: HarpoTrackingPh.h:97
HarpoDccMap::print
virtual void print()
Definition: HarpoDccMap.cxx:124
HarpoEventHeader
FullEvent Header not scecific to the detectors The class is ....
Definition: HarpoEventHeader.h:27
HarpoTrackingPh::print
void print()
Ovreloaded medod whic do all job.
Definition: HarpoTrackingPh.cxx:41
HarpoTrackingPh::GetMapEdges
Int_t GetMapEdges(TClonesArray *clArray, Int_t plane, Int_t &iMin, Int_t &iMax, Int_t &jMin, Int_t &jMax)
Definition: HarpoTrackingPh.cxx:307
HarpoRecoClusters::RemoveAllClusterTrack
void RemoveAllClusterTrack()
Definition: HarpoRecoEvent.h:74
HarpoTrackingPh::GetQtrack
Double_t GetQtrack(Int_t itr, Int_t plane)
Definition: HarpoTrackingPh.cxx:865
HarpoEventHeader::print
virtual void print()
Definition: HarpoEventHeader.cxx:21
HarpoTrackingPh::fNclTrack
Int_t fNclTrack[NTRACK]
Definition: HarpoTrackingPh.h:83
HarpoRecoKalmanTracks
HarpoRecoTracks object, Obtained with Kalman filter.
Definition: HarpoRecoEvent.h:223
HarpoTrackingPh::used
Int_t used[4000]
Definition: HarpoTrackingPh.h:107
HarpoRecoClusters
Cluster object, containing position, charge and quality information.
Definition: HarpoRecoEvent.h:28
HarpoAnalyse::OpenHistFile
TFile * OpenHistFile(const char *ananame)
Definition: HarpoAnalyse.cxx:150
type
static int type
Definition: HarpoConfig.cxx:106
HarpoDccMap
unpacked dcc data The class contains the data map for DCC or Feminos The data is stored as a 2D TMatr...
Definition: HarpoDccMap.h:29
HarpoTrackingPh::GetQtracks
Double_t GetQtracks(Int_t plane)
Definition: HarpoTrackingPh.cxx:903
HarpoRecoClusters::GetIdClusterTrack
Int_t GetIdClusterTrack()
Definition: HarpoRecoEvent.h:45
HarpoTrackingPh::fTrack
Int_t fTrack[4000]
Definition: HarpoTrackingPh.h:109
HarpoTrackingPh::InitPlane
Int_t InitPlane(TClonesArray *clArray, Int_t plane)
Definition: HarpoTrackingPh.cxx:210
HarpoTrackingPh::FindNextClosest
TArrayI * FindNextClosest(TMatrixD X, TMatrixD C, Double_t angle, Int_t Ntr, TClonesArray *clArray, TArrayI *arr, Int_t ncl, Int_t plane, Int_t color, Int_t fill, Int_t smooth, Int_t q)
Definition: HarpoTrackingPh.cxx:992
HarpoRecoClusters::SetXfit
void SetXfit(Int_t i, Double_t val)
Definition: HarpoRecoEvent.h:70
HarpoTrackingPh::reused
Int_t reused[4000]
Definition: HarpoTrackingPh.h:108
NADC
#define NADC
Definition: HarpoDccMap.h:18
KALMAN
#define KALMAN
Definition: HarpoRecoEvent.h:26
HarpoTrackingPh::FindNext
TArrayI * FindNext(TMatrixD X, TMatrixD C, Double_t angle, Int_t Ntr, TClonesArray *clArray, TArrayI *arr, Int_t ncl, Int_t plane, Int_t color, Int_t fill, Int_t smooth, Int_t q)
Definition: HarpoTrackingPh.cxx:985
HarpoTrackingPh::fStartPointX
Double_t fStartPointX[NTRACK][10]
Definition: HarpoTrackingPh.h:87
TCLUSTER
#define TCLUSTER
Definition: HarpoRecoEvent.h:23
gHarpoDebug
Long64_t gHarpoDebug
Definition: HarpoDebug.cxx:9
HarpoTrackingPh::hThetaDist
TH2F * hThetaDist[10]
Definition: HarpoTrackingPh.h:117
HarpoTrackingPh::fGraphs
TGraph * fGraphs[2][NTRACK]
Definition: HarpoTrackingPh.h:84
HarpoTrackingPh::fEndIndex
Int_t fEndIndex[NTRACK]
Definition: HarpoTrackingPh.h:98
HarpoTrackingPh::fQcommon
Double_t fQcommon[NTRACK][NTRACK]
Definition: HarpoTrackingPh.h:100
HarpoEvent::GetHeader
HarpoEventHeader * GetHeader()
Definition: HarpoEvent.cxx:80
HarpoTrackingPh::fMaxSlopeType
Double_t fMaxSlopeType
Definition: HarpoTrackingPh.h:136
HarpoAnalyse::fEvt
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
HarpoTrackingPh::NTcl
Int_t NTcl[NADC]
Definition: HarpoTrackingPh.h:80
HarpoTrackingPh::Tcl
Int_t Tcl[NADC][20]
Definition: HarpoTrackingPh.h:79
HarpoAnalyse::nEvents
ULong_t nEvents
Definition: HarpoAnalyse.h:75
HarpoTrackingPh::AddTrack
Int_t AddTrack(TClonesArray *clArray, Int_t plane)
Definition: HarpoTrackingPh.cxx:414
res
unsigned int res
Definition: Pmm2Status.h:428
HarpoTrackingPh::fStartDirZ
Double_t fStartDirZ[NTRACK][10]
Definition: HarpoTrackingPh.h:91
gkNDetectors
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
CCLUSTER
#define CCLUSTER
Definition: HarpoRecoEvent.h:22
HarpoEvent::GetDccMap
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
HarpoTrackingPh
A virtual class which define intrafece between HARPO Reader and Event Analysis code.
Definition: HarpoTrackingPh.h:26
HarpoTrackingPh::fMapTmp
ULong64_t fMapTmp[NCHAN][NADC]
Definition: HarpoTrackingPh.h:103
HarpoRecoClusters::SetZfit
void SetZfit(Int_t i, Double_t val)
Definition: HarpoRecoEvent.h:71
HarpoTrackingPh::Ccl
Int_t Ccl[NCHAN][20]
Definition: HarpoTrackingPh.h:77
HarpoRecoEvent::GetClustersArray
TClonesArray * GetClustersArray()
Definition: HarpoRecoEvent.h:594
HarpoTrackingPh::GetResolution
Double_t GetResolution(HarpoRecoClusters *cl)
Definition: HarpoTrackingPh.cxx:1303
HarpoRecoEvent::SetTrackType
void SetTrackType(Int_t val)
Definition: HarpoRecoEvent.h:663
HarpoTrackingPh::hNcl
TH1F * hNcl
Redefine empty default.
Definition: HarpoTrackingPh.h:39
HarpoRecoClusters::CheckIdClusterTrack
Bool_t CheckIdClusterTrack(Int_t val)
Definition: HarpoRecoEvent.h:79
Generated by   doxygen 1.8.6