HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoHoughTracking.cxx
Go to the documentation of this file.
1 //
2 // File HarpoHoughTracking.cxx
3 //
12 #include "HarpoHoughTracking.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 #include "MakeNiceHisto.h"
21 
22 #include "TFile.h"
23 #include "TStyle.h"
24 #include "TCanvas.h"
25 #include "TRootEmbeddedCanvas.h"
26 #include "TLatex.h"
27 #include "TGraph.h"
28 #include "TF1.h"
29 #include "TLinearFitter.h"
30 #include "TLatex.h"
31 #include "TGLabel.h"
32 #include "TMath.h"
33 
34 #include <cstdlib>
35 #include <cstring>
36 #include <cassert>
37 #include <fstream>
38 #include <iostream>
39 using namespace std;
40 
41 ClassImp(HarpoHoughTracking);
42 
43 const Double_t HarpoHoughTracking::kfWidthRho = 200;
44 const Double_t HarpoHoughTracking::kfWidthTheta = 4;
45 const Double_t HarpoHoughTracking::kfFacHough = 4;
46 
47 void HarpoHoughTracking::print()
48 {
49  unsigned int nd; // number of detectors
50  HarpoEventHeader *hdr = fEvt->GetHeader();
51  assert(hdr != NULL);
52  hdr->print();
53 
54  for (nd = 0; nd < gkNDetectors; nd++) {
55  // if (fEvt->isdataExist(nd)) {
56  HarpoDccMap *plane = fEvt->GetDccMap(nd);
57  if (plane != NULL )plane->print();
58  }
59 }
60 
61 void HarpoHoughTracking::process()
62 {
63  nEvents++;
64  if(gHarpoDebug>0 || nEvents%1000 == 0)
65  Info("process","Process event %ld",nEvents);
66  // if(fChooseQmin) SetConfig();
67  if(fEvt->GetRecoEvent()->GetTracksArray()->GetEntries()>0){
68  Info("process","Reset existing tracking %ld",nEvents);
69  fEvt->GetRecoEvent()->ResetTracksArray();
70  TClonesArray* clArray = fEvt->GetRecoEvent()->GetClustersArray();
71  Int_t nCl = clArray->GetEntries();
72  for(Int_t i=0;i<nCl;i++){
73  HarpoRecoClusters* cl1 = (HarpoRecoClusters*)clArray->At(i);
74  cl1->RemoveAllClusterTrack();
75  }
76  }
77  fEvt->GetRecoEvent()->SetTrackType(HOUGH);
78  if(fDisplay){
79  hSigmaRhoEvent->Reset();
80  hSigmaThetaEvent->Reset();
81  }
82  InitTrack();
83  // if(fCanvas){
84  // // fCanvas->Divide(2);
85  // fCanvas->GetPad(1)->Delete();
86  // fCanvas->GetPad(2)->Delete();
87  // fCanvas->Divide(2);
88  // fCanvas->GetPad(1)->Divide(5,5);
89  // fCanvas->GetPad(2)->Divide(5,5);
90  // }
91  for(Int_t ndet=0;ndet<2;ndet++) { // Detectors 0 and 1 = X and Y planes
92  if (!gHDetSet->isExist(ndet)) continue;
93  //MAKE TRACK
94  FindTrack(ndet);
95  }
96  // if(fCanvas){
97  // fCanvas->GetPad(1)->cd(25);
98  // hSigmaRho->Draw();
99  // hSigmaTheta->Draw("same");
100  // fCanvas->Update();
101  // }
102 }
103 
104 
105 
106 
107 
108 
109 void HarpoHoughTracking::FindPairNew(Int_t plane, Int_t trackId, Double_t &muRhoOut, Double_t &muThetaOut, Double_t &sigRhoOut, Double_t &sigThetaOut)
110 {
111  // if(gHarpoDebug>1)
112  // Info("FindPairNew","************ %f %d %d %d %d ******************",cutHough, cutnclmin, cutnclmax, cutqmin, plane);
113  //Info("FindPairNew","******************************");
114  double pi=TMath::Pi();
115 
116  Double_t rapt0 = 0;
117  Double_t rapr0 = 0;
118  Int_t nloop=0;
119  Int_t goodbinning=1;
120 
121  Double_t xCenter = 0.;
122  Double_t yCenter = 0.;
123  Double_t xWidth = kfWidthRho;
124  Double_t yWidth = kfWidthTheta;
125 
126  Int_t maxLoops = 6;
127  while(goodbinning==1 && nloop < maxLoops){
128  //for(Int_t k = 0; k<5; k++){
129  //if(nloop == maxLoops - 1)
130  if(gHarpoDebug>1)
131  Info("FindPairNew","loop == %d",nloop);
132 
133  TH2F* h = new TH2F("h",";#rho [mm];#theta",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
134  TH2F* hPair = new TH2F("hPair",";#rho [mm];#theta",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
135  TH2F* hPairW = new TH2F("hPairW",";#rho [mm];#theta",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
136  TH2F* hPairRho = new TH2F("hPairRho","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
137  TH2F* hPairRho2 = new TH2F("hPairRho2","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
138  TH2F* hPairRho3 = new TH2F("hPairRho3","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
139  TH2F* hPairRho4 = new TH2F("hPairRho4","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
140  TH2F* hPairTheta = new TH2F("hPairTheta","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
141  TH2F* hPairTheta2 = new TH2F("hPairTheta2","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
142  TH2F* hPairTheta3 = new TH2F("hPairTheta3","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
143  TH2F* hPairTheta4 = new TH2F("hPairTheta4","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
144  TH2F* hPairRhoTheta = new TH2F("hPairRhoTheta","",fNbinr,xCenter - xWidth, xCenter + xWidth,fNbint,yCenter - yWidth, yCenter + yWidth);
145 
146 
147  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
148  TClonesArray* clArray = reco->GetClustersArray();
149  Int_t nCl = clArray->GetEntries();
150  for(Int_t i=0;i<nCl;i++){
151  HarpoRecoClusters* cl1 = (HarpoRecoClusters*)clArray->At(i);
152  // Int_t ind1 = cl1->GetIndex();
153  // Double_t mean1 = cl1->GetMean();
154  Int_t type1 = cl1->GetType();
155  Double_t q1 = cl1->GetQ();
156 
157  //Info("FindPairNew","%d %d %d %d %d",cl1->GetPlane(),iclusttrack[i],cl1->GetQ(),cl1->GetWidth(),cl1->GetType());
158  if(cl1->GetPlane()!=plane) continue;
159  if(iclusttrack[i]>=0) continue;
160  if(cl1->GetQ()<=fQmin[type1]) continue;
161  if(cl1->GetWidth()>=fWmax[type1]) continue;
162  if(cl1->GetWidth()<=fWmin[type1]) continue;
163  //Info("FindPairNew","Cluster accepted");
164 
165  Double_t t1 = 0, c1 = 0;
166  t1 = cl1->GetZ() - NADC/2;
167  c1 = cl1->GetX() - NCHAN/2;
168  // if(type1==TCLUSTER){
169  // t1 = ind1 - NADC/2;
170  // c1 = mean1 - NCHAN/2;
171  // }
172  // if(type1==CCLUSTER){
173  // t1 = mean1 - NADC/2;
174  // c1 = ind1 - NCHAN/2;
175  // }
176  for(Int_t j=i+1;j<nCl;j++){
177  HarpoRecoClusters* cl2 = (HarpoRecoClusters*)clArray->At(j);
178  // Int_t ind2 = cl2->GetIndex();
179  // Double_t mean2 = cl2->GetMean();
180  Int_t type2 = cl2->GetType();
181  Double_t q2 = cl2->GetQ();
182  if(cl2->GetPlane()!=plane) continue;
183  if(iclusttrack[j]>=0) continue;
184  if(cl2->GetQ()<=fQmin[type2]) continue;
185  if(cl2->GetWidth()>=fWmax[type2]) continue;
186  if(cl2->GetWidth()<=fWmin[type2]) continue;
187 
188  Double_t t2 = 0, c2 = 0;
189  t2 = cl2->GetZ() - NADC/2;
190  c2 = cl2->GetX() - NCHAN/2;
191  // if(type2==TCLUSTER){
192  // t2=ind2 - NADC/2;
193  // c2=mean2 - NCHAN/2;
194  // }
195  // if(type2==CCLUSTER){
196  // t2=mean2 - NADC/2;
197  // c2=ind2 - NCHAN/2;
198  // }
199  Double_t dt = t2 - t1;
200  Double_t dc = c2 - c1;
201  Double_t dist = sqrt(dt*dt + dc*dc);
202 
203  if (dist<=fDistMin) continue;
204  Double_t theta = pi/2;
205  if(dt) theta = -TMath::ATan(dc/dt);
206  Double_t cross = c1*t2 - c2*t1;
207  Double_t rho = TMath::Abs(cross/dist);
208  if((t2-t1)*cross<0) rho = -rho;
209  if(type1 == CCLUSTER && type2 == CCLUSTER && TMath::Abs(TMath::Cos(theta)) < 0.5) continue;
210  if(type1 == TCLUSTER && type2 == TCLUSTER && TMath::Abs(TMath::Sin(theta)) < 0.5) continue;
211  Double_t weight = 1;
212  if(q1>0 && q2>0) weight = TMath::Sqrt(q1*q2);
213  //Info("FindPairNew","rho=%g, theta=%g, weight=%g, %d %d",rho,theta,weight,q1,q2);
214  hPair->Fill(rho, theta);
215  hPairW->Fill(rho, theta,weight);
216  hPairRho->Fill(rho, theta, rho*weight);
217  hPairRho2->Fill(rho, theta, rho*rho*weight);
218  hPairTheta->Fill(rho, theta, theta*weight);
219  hPairTheta2->Fill(rho, theta, theta*theta*weight);
220  hPairRhoTheta->Fill(rho, theta, rho*theta*weight);
221 
222  // Wrap angle around (modulo pi)
223  rho = -rho;
224  if(theta>0) theta = theta - TMath::Pi();
225  else theta = theta + TMath::Pi();
226  hPair->Fill(rho, theta);
227  hPairW->Fill(rho, theta,weight);
228  hPairRho->Fill(rho, theta, rho*weight);
229  hPairRho2->Fill(rho, theta, rho*rho*weight);
230  hPairRho3->Fill(rho, theta, rho*rho*rho*weight);
231  hPairRho4->Fill(rho, theta, rho*rho*rho*rho*weight);
232  hPairTheta->Fill(rho, theta, theta*weight);
233  hPairTheta2->Fill(rho, theta, theta*theta*weight);
234  hPairTheta3->Fill(rho, theta, theta*theta*theta*weight);
235  hPairTheta4->Fill(rho, theta, theta*theta*theta*theta*weight);
236  hPairRhoTheta->Fill(rho, theta, rho*theta*weight);
237  }
238  }
239 
240 
241  //Info("FindPairNew","N = %.0f",hPair->GetEntries());
242  if(hPair->GetEntries()<10) {
243  h->Delete();
244  hPair->Delete();
245  hPairW->Delete();
246  hPairRho->Delete();
247  hPairRho2->Delete();
248  hPairRho3->Delete();
249  hPairRho4->Delete();
250  hPairTheta->Delete();
251  hPairTheta2->Delete();
252  hPairTheta3->Delete();
253  hPairTheta4->Delete();
254  hPairRhoTheta->Delete();
255  break;
256  }
257 
258  Double_t sRmax = 0;//, sTmax = 0;
259  // Calculate Mean and RMS in histograms
260  for (Int_t ir=1;ir<=fNbinr;ir++){
261  for (Int_t it=1;it<=fNbint;it++){
262  Int_t n = hPair->GetBinContent(ir,it);
263  Double_t weight = hPairW->GetBinContent(ir,it);
264  if(n<=0) continue;
265  if(weight<=0) continue;
266  Double_t muRhoTmp = hPairRho->GetBinContent(ir,it)/weight;
267  Double_t m3RhoTmp = hPairRho3->GetBinContent(ir,it)/weight;
268  Double_t m4RhoTmp = hPairRho4->GetBinContent(ir,it)/weight;
269  Double_t sigmaRhoTmp = hPairRho2->GetBinContent(ir,it)/weight - muRhoTmp*muRhoTmp;
270  m4RhoTmp = m4RhoTmp - 4*muRhoTmp*m3RhoTmp + 6*muRhoTmp*muRhoTmp*sigmaRhoTmp - 3*muRhoTmp*muRhoTmp*muRhoTmp*muRhoTmp;
271  Double_t muThetaTmp = hPairTheta->GetBinContent(ir,it)/weight;
272  Double_t m3ThetaTmp = hPairTheta3->GetBinContent(ir,it)/weight;
273  Double_t m4ThetaTmp = hPairTheta4->GetBinContent(ir,it)/weight;
274  // Double_t sigmaThetaTmp = hPairTheta2->GetBinContent(ir,it)/weight - muThetaTmp*muThetaTmp;
275  m4ThetaTmp = m4ThetaTmp - 4*muThetaTmp*m3ThetaTmp + 6*muThetaTmp*muThetaTmp*sigmaRhoTmp - 3*muThetaTmp*muThetaTmp*muThetaTmp*muThetaTmp;
276  if(sigmaRhoTmp > sRmax) sRmax = sigmaRhoTmp;
277  // if(sigmaRhoTmp>1e-5 && sigmaThetaTmp>1e-5 && (n>0)){
278  // if(sigmaRhoTmp>0 && sigmaThetaTmp>0 && n>2){
279  // Double_t vRho = (m4RhoTmp - (n-3.)*sigmaRhoTmp*sigmaRhoTmp/(n-1))/n;
280  // Double_t vTheta = (m4ThetaTmp - (n-3.)*sigmaThetaTmp*sigmaThetaTmp/(n-1))/n;
281  // h->SetBinContent(ir,it,(n-1.)/TMath::Sqrt(sigmaRhoTmp*sigmaThetaTmp));
282  // //h->SetBinContent(ir,it,vRho/sigmaRhoTmp + vTheta/sigmaThetaTmp);
283  // }else h->SetBinContent(ir,it,0);
284  // if(hPair->GetMaximum()<=500)
285  h->SetBinContent(ir,it,n);
286  }
287  }
288 
289 
290 
291 
292  // Convert to 1D histograms around the maximum bin
293  Int_t ix = 0, iy = 0, iz = 0;
294  h->GetYaxis()->SetRangeUser(-TMath::Pi()/2,TMath::Pi()/2);
295  h->GetMaximumBin(ix, iy, iz);
296  Double_t rhoMax = hPair->GetXaxis()->GetBinLowEdge(ix);
297  Double_t thetaMax = hPair->GetYaxis()->GetBinLowEdge(iy);
298 
299 
300  if(gHarpoDebug>1)
301  Info("","%g %g %g",hPairRho2->GetBinContent(ix,iy),hPairTheta2->GetBinContent(ix,iy),hPair->GetBinContent(ix,iy));
302 
303  Double_t muRhoPeak = hPairRho->GetBinContent(ix,iy)/hPairW->GetBinContent(ix,iy);
304  Double_t sigmaRhoPeak = TMath::Sqrt(hPairRho2->GetBinContent(ix,iy)/hPairW->GetBinContent(ix,iy) - muRhoPeak*muRhoPeak)/hPair->GetXaxis()->GetBinWidth(1);
305  Double_t muThetaPeak = hPairTheta->GetBinContent(ix,iy)/hPairW->GetBinContent(ix,iy);
306  Double_t sigmaThetaPeak = TMath::Sqrt(hPairTheta2->GetBinContent(ix,iy)/hPairW->GetBinContent(ix,iy) - muThetaPeak*muThetaPeak)/hPair->GetYaxis()->GetBinWidth(1);
307 
308 
309  //Info("FindPairNew","Draw %d %d %d | %g %g %g | %.0f",plane,trackId,nloop,sigmaRhoPeak,sigmaThetaPeak,hPair->GetBinContent(ix,iy),hPair->GetEntries());
310  // if(fCanvas && trackId<4 && nloop<5){
311  // Info("FindPairNew","Draw %d %d %d | %g %g %d",plane,trackId,nloop,sigmaRhoPeak,sigmaThetaPeak,hPair->GetBinContent(ix,iy));
312  // MakeNiceHisto(h,fCanvas->GetPad(plane+1)->GetPad(5*(trackId) + nloop + 1),"col");
313  // // if(fCanvas->GetPad(plane+1)->GetPad(5*trackId+1) + nloop + 1))
314  // // fCanvas->GetPad(plane+1)->GetPad(5*trackId+1) + nloop + 1)->SetLogz();
315  // TLatex* latex = new TLatex();
316  // latex->SetTextFont(132);
317  // latex->SetTextAlign(12);
318  // fCanvas->GetPad(plane+1)->cd(5*trackId + nloop + 1);
319  // latex->DrawLatex(xCenter + 0.5*xWidth,yCenter+0.5*yWidth,Form("TR%d",trackId));
320  // }
321  if(gHarpoDebug>1)
322  Info("FindPairNew","###################### %d %d %d (%d)",plane,trackId,nloop,fDisplay);
323 
324  if(fDisplay && trackId<5 && nloop<5){
325  Info("FindPairNew","%d %d %d => %g %g %g",plane,trackId,nloop,sigmaRhoPeak,sigmaThetaPeak,hPair->GetBinContent(ix,iy));
326  fHist[plane][trackId][nloop] = (TH2F*)h->Clone(Form("h%d_%d_%d",plane,trackId,nloop));
327  if(nloop>fNloopMax) fNloopMax = nloop;
328  if(trackId>fTrackIdMax) fTrackIdMax = trackId;
329  fSigmaRhoPeak[plane][trackId][nloop] = sigmaRhoPeak;
330  fSigmaThetaPeak[plane][trackId][nloop] = sigmaThetaPeak;
331  fNpairsPeak[plane][trackId][nloop] = hPair->GetBinContent(ix,iy);
332  }
333 
334  hSigmaRho->Fill(sigmaRhoPeak);
335  hSigmaTheta->Fill(sigmaThetaPeak);
336  hSigmaRhoEvent->Fill(sigmaRhoPeak);
337  hSigmaThetaEvent->Fill(sigmaThetaPeak);
338  if(gHarpoDebug>1)
339  Info("FindPairNew","sigmaRhoPeak = %g, sigmaThetaPeak = %g, hPair = %g",sigmaRhoPeak,sigmaThetaPeak,hPair->GetXaxis()->GetBinWidth(1));
340  if((sigmaRhoPeak<fSigmaMaxRebin || sigmaThetaPeak<fSigmaMaxRebin) && nloop < maxLoops - 1 && hPair->GetBinContent(ix,iy)>fMinPairs){
341  // Binning is not fine enough => refine binning around max
342  goodbinning = 1;
343  xCenter = rhoMax;
344  yCenter = thetaMax;
345 
346  xWidth /= kfFacHough;
347  yWidth /= kfFacHough;
348  nloop++;
349  h->Delete();
350  hPair->Delete();
351  hPairW->Delete();
352  hPairRho->Delete();
353  hPairRho2->Delete();
354  hPairRho3->Delete();
355  hPairRho4->Delete();
356  hPairTheta->Delete();
357  hPairTheta2->Delete();
358  hPairTheta3->Delete();
359  hPairTheta4->Delete();
360  hPairRhoTheta->Delete();
361  continue;
362  }
363 
364 
365 
366 
367 
368 
369 
370 
371 
372 
373 
374 
375 
376 
377  TH1F* h1D = new TH1F("h1D","",fNbinr,0,fNbinr);
378  TH1F* h1DRho = new TH1F("h1DRho","",fNbinr,0,fNbinr);
379  TH1F* h1DRho2 = new TH1F("h1DRho2","",fNbinr,0,fNbinr);
380  TH1F* h1DTheta = new TH1F("h1DTheta","",fNbinr,0,fNbinr);
381  TH1F* h1DTheta2 = new TH1F("h1DTheta2","",fNbinr,0,fNbinr);
382  TH1F* h1DRhoTheta = new TH1F("h1DRhoTheta","",fNbinr,0,fNbinr);
383  for (Int_t ir=1;ir<=fNbinr;ir++){
384  for (Int_t it=1;it<=fNbint;it++){
385  Int_t idx=TMath::Abs(ir-ix);
386  Int_t idy=TMath::Abs(it-iy);
387  Double_t idd = idx + idy + 0.5;
388  //if(idd == 1.5) Info("FindPairNew","%d %d %g",ir,it,hPair->GetBinContent(ir,it));
389  h1D->Fill(idd,hPairW->GetBinContent(ir,it));
390  h1DRho->Fill(idd,hPairRho->GetBinContent(ir,it));
391  h1DRho2->Fill(idd,hPairRho2->GetBinContent(ir,it));
392  h1DTheta->Fill(idd,hPairTheta->GetBinContent(ir,it));
393  h1DTheta2->Fill(idd,hPairTheta2->GetBinContent(ir,it));
394  h1DRhoTheta->Fill(idd,hPairRhoTheta->GetBinContent(ir,it));
395 
396  // Double_t muRhoTmp = hPairRho->GetBinContent(ir,it)/hPair->GetBinContent(ir,it);
397  // Double_t sigmaRhoTmp = hPairRho2->GetBinContent(ir,it)/hPair->GetBinContent(ir,it) - muRhoTmp*muRhoTmp;
398  // Double_t muThetaTmp = hPairTheta->GetBinContent(ir,it)/hPair->GetBinContent(ir,it);
399  // Double_t sigmaThetaTmp = hPairTheta2->GetBinContent(ir,it)/hPair->GetBinContent(ir,it) - muThetaTmp*muThetaTmp;
400  // if(sigmaRhoTmp>0 && sigmaThetaTmp>0 && hPair->GetBinContent(ir,it)>0)
401  // hPairRho2->SetBinContent(ir,it,(hPair->GetBinContent(ir,it)-1.)/TMath::Sqrt(sigmaRhoTmp*sigmaThetaTmp));
402  // else hPairRho2->SetBinContent(ir,it,0);
403  }
404  }
405 
406 
407 
408 
409 
410  Int_t ioMax = h1D->GetXaxis()->GetNbins();
411  for(Int_t io=2;io<=ioMax;io++){
412  Int_t n = h1D->Integral(1,io);
413  Double_t muRho=h1DRho->Integral(1,io)/n;
414  Double_t sigRho=h1DRho2->Integral(1,io)/n - muRho*muRho;
415 
416  if(sigRho>0) sigRho = TMath::Sqrt(sigRho);
417  else sigRho = 0;
418 
419  Double_t muTheta=h1DTheta->Integral(1,io)/n;
420 
421  Double_t sigTheta=h1DTheta2->Integral(1,io)/n - muTheta*muTheta;
422  if(sigTheta>0) sigTheta = sqrt(sigTheta);
423  else sigTheta = 0;
424 
425  // Double_t covRhoTheta = h1DRhoTheta->Integral(1,io)/n - muRho * muTheta;
426  // Double_t corRhoTheta = 0;
427  // if(sigRho>0&&sigTheta>0) corRhoTheta = covRhoTheta / (sigRho * sigTheta);
428  // Double_t rap=0;
429  Double_t rapt=0;
430  Double_t rapr=0;
431  if(sigRho>0) rapr = (io-1) * hPair->GetXaxis()->GetBinWidth(1) / sigRho;
432  if(sigTheta>0) rapt = (io-1) * hPair->GetYaxis()->GetBinWidth(1) / sigTheta;
433  // if(rapt>0) rap = rapr/rapt;
434 
435  if(gHarpoDebug>1)
436  Info("FindPairNew","io = %d, rapr = %g, rapr0 = %g, rapt = %g, rapt0 = %g", io,rapr,rapr0,rapt,rapt0);
437  if(io==ioMax || (((rapr>fCutoffHough && rapr0<fCutoffHough)||(rapt>fCutoffHough && rapt0<fCutoffHough))&&io>2)){
438  goodbinning = 0;
439  muRhoOut = muRho;
440  sigRhoOut = sigRho;
441  muThetaOut = muTheta;
442  sigThetaOut = sigTheta;
443 
444  if(gHarpoDebug>1 || fDisplay)
445  Info("FindPairNew","theta = %.3g +- %.3g rho = %.3g +- %.3g",muThetaOut, sigThetaOut, muRhoOut, sigRhoOut);
446  //break;
447  h->Delete();
448  hPair->Delete();
449  hPairW->Delete();
450  hPairRho->Delete();
451  hPairRho2->Delete();
452  hPairRho3->Delete();
453  hPairRho4->Delete();
454  hPairTheta->Delete();
455  hPairTheta2->Delete();
456  hPairTheta3->Delete();
457  hPairTheta4->Delete();
458  hPairRhoTheta->Delete();
459  h1D->Delete();
460  h1DRho->Delete();
461  h1DRho2->Delete();
462  h1DTheta->Delete();
463  h1DTheta2->Delete();
464  h1DRhoTheta->Delete();
465  return;
466  }
467  rapr0 = rapr;
468  rapt0 = rapt;
469 
470 
471  }
472  h->Delete();
473  hPair->Delete();
474  hPairW->Delete();
475  hPairRho->Delete();
476  hPairRho2->Delete();
477  hPairRho3->Delete();
478  hPairRho4->Delete();
479  hPairTheta->Delete();
480  hPairTheta2->Delete();
481  hPairTheta3->Delete();
482  hPairTheta4->Delete();
483  hPairRhoTheta->Delete();
484  h1D->Delete();
485  h1DRho->Delete();
486  h1DRho2->Delete();
487  h1DTheta->Delete();
488  h1DTheta2->Delete();
489  h1DRhoTheta->Delete();
490  // Info("FindPairNew","** loop = %d (%d)",nloop,goodbinning);
491  }
492 
493 }
494 
495 
496 
497 
498 
499 
500 
501 
502 
503 
504 
505 
506 void HarpoHoughTracking::FindTrack(Int_t plane)
507 {
508  if(gHarpoDebug>1) Info("FindTrack","Plane %d",plane);
509  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
510  TClonesArray* clArray = reco->GetClustersArray();
511  Int_t nCl = clArray->GetEntries();
512  // Double_t stheta,ctheta,sigtrans;
513  // Int_t rlambda;
514  // Double_t dist;
515  // Double_t x0,y0,x,y;
516  Int_t mtrack = 0;
517  Int_t nClLeftOld = -1;
518 
519  Int_t trackId = 0;
520  do {
521  Int_t Ncluster=0;
522  Double_t rho = -1000, theta = -1000,sigRho = -1, sigTheta = -1;
523  FindPairNew(plane, trackId,rho,theta,sigRho,sigTheta);
524  if(theta == -1000) break;
525  if(sigTheta<1e-3) sigTheta = 0.001;
526  if(sigRho<1e-3) sigRho = 0.001;
527  Double_t stheta = cos(theta);
528  Double_t ctheta = sin(theta);
529 
530  // Int_t j;
531  Int_t ic;
532  Int_t typeTrack = -1;//ttz = -1;
533  Int_t nClOnly = 0, nClLeft = 0;
534  for (ic=0;ic<nCl;ic++){
535  HarpoRecoClusters* cl = (HarpoRecoClusters*)clArray->At(ic);
536  Int_t type = cl->GetType();
537  if(type == CCLUSTER && TMath::Abs(ctheta) < 0.5) continue;
538  if(type == TCLUSTER && TMath::Abs(stheta) < 0.5) continue;
539  Int_t planeCl = cl->GetPlane();
540  // Double_t q = cl->GetQ();
541  if(planeCl!=plane) continue;
542  //if(iclusttrack[ic]>=0) continue;
543  if(cl->GetQ()<=fQmin[type]) continue;
544  if(cl->GetWidth()>=fWmax[type]) continue;
545  if(cl->GetWidth()<=fWmin[type]) continue;
546  if(gHarpoDebug>1)
547  Info("FindTrack","cluster %d, Plane= %d", ic, plane);
548 
549  Double_t x0 = cl->GetZ();
550  Double_t y0 = cl->GetX();
551 
552 
553  Double_t rlambda = stheta*(x0 - NADC/2) - ctheta*(y0 - NCHAN/2) ;
554  Double_t x = rho*ctheta + rlambda*stheta + NADC/2;
555  Double_t y = rho*stheta - rlambda*ctheta + NCHAN/2;
556  //calcul de distance d'un cluster à une trace(theta et rho):
557  Double_t dist = (x - x0)*ctheta + (y - y0)*stheta ;
558  Double_t sigtrans = sigRho;//sqrt(sigRho*sigRho+rlambda*sigTheta*rlambda*sigTheta);
559 
560  if(cl->GetNtr()<1) nClLeft++;
561  if(TMath::Abs(dist)< RouteHoughLoose*sigtrans){
562  Ncluster++;
563  if(cl->GetNtr()<1) nClOnly++;
564  }
565  if(TMath::Abs(dist)< RouteHoughTight*sigtrans)
566  iclusttrack[ic] = 1;
567 
568  }
569 
570  if(gHarpoDebug>0 || fDisplay)
571  Info("FindTrack","Ncluster = %d/%d/%d, Plane= %d",nClOnly, Ncluster,nClLeft, plane);
572 
573  if ( nClLeft == nClLeftOld) break;
574  nClLeftOld = nClLeft;
575  // if(nClOnly<0.4*Ncluster || nClOnly<5) continue;
576  if(nClOnly<5) continue;
577 
578  if( Ncluster<=fNclMin)
579  continue;
580  trackId++;
581 
582  Ncluster = 0;
583 
584  TGraph* gTrack = new TGraph();
585  Int_t nDraw = 10;
586  Double_t lambdaMin = 1e8, lambdaMax = -1e8;
587  Double_t xStart = 0, zStart = 0, xEnd = 0, zEnd = 0;
588  for (ic=0;ic<nCl;ic++){
589  HarpoRecoClusters* cl = (HarpoRecoClusters*)clArray->At(ic);
590  // Int_t ind = cl->GetIndex();
591  // Double_t mean = cl->GetMean();
592  // Int_t planeCl = cl->GetPlane();
593  Double_t q = cl->GetQ();
594  Int_t type = cl->GetType();
595  // Double_t qMax = cl-> GetQmax();
596  // Double_t sigma = cl->GetSig();
597  // Int_t width= cl->GetWidth();
598  if(cl->GetPlane()!=plane) continue;
599  if(cl->GetQ()<=fQmin[type]) continue;
600  if(cl->GetWidth()>=fWmax[type]) continue;
601  if(cl->GetWidth()<=fWmin[type]) continue;
602 
603  Double_t x0 = cl->GetZ();
604  Double_t y0 = cl->GetX();
605 
606  Double_t rlambda = stheta*(x0 - NADC/2) - ctheta*(y0 - NCHAN/2) ;
607  Double_t x = rho*ctheta + rlambda*stheta + NADC/2;
608  Double_t y = rho*stheta - rlambda*ctheta + NCHAN/2;
609  //calcul de distance d'un cluster à une trace(theta et rho):
610  Double_t dist = (x - x0)*ctheta + (y - y0)*stheta ;
611  Double_t sigtrans = sigRho;//sqrt(sigRho*sigRho+rlambda*sigTheta*rlambda*sigTheta);
612 
613  nDraw++;
614  if(TMath::Abs(dist)< RouteHoughTight*sigtrans)
615  iclusttrack[ic] = 1;
616  if(TMath::Abs(dist)< RouteHoughLoose*sigtrans){
617  if(type==TCLUSTER) {
618  // j=cl->GetIndex();
619  Qtrack=Qtrack+ q;
620  }
621  if(type==BLOCCLUSTER)
622  Qtrack=Qtrack+ q;
623 
624  if(1){//cl->GetNtr()<1){
625  cl->AddIdClusterTrack(mtrack);
626  gTrack->SetPoint(gTrack->GetN(),x0,y0);
627  Ncluster++;
628  if(rlambda < lambdaMin){
629  lambdaMin = rlambda;
630  xStart = y0;
631  zStart = x0;
632  }
633  if(rlambda > lambdaMax){
634  lambdaMax = rlambda;
635  xEnd = y0;
636  zEnd = x0;
637  }
638  }
639  }
640  }
641  Double_t slope = -1000, x0 = -1000;
642  if(stheta){
643  slope=-ctheta/stheta;
644  x0 = NCHAN/2 + rho/stheta-slope*NADC/2;
645  }
646 
647  // if(gTrack->GetN()>10){
648  // TLinearFitter* linfit = new TLinearFitter(1,"pol1","D");
649  // linfit->AssignData(gTrack->GetN(), 1, gTrack->GetX(), gTrack->GetY());
650  // linfit->EvalRobust();
651  // x0 = linfit->GetParameter(0);
652  // slope = linfit->GetParameter(1);
653  // // chi2 = linfit->GetChisquare();
654  // linfit->Delete();
655  // }
656 
657  // Double_t pi =TMath::Pi();
658  Double_t AngleTrack = 0.0;
659  // AngleTrack = -Mtheta;
660  AngleTrack = slope;//-theta;
661  if(gHarpoDebug>0)
662  Info("FindTrack","%g / %g = %g",ctheta,stheta,AngleTrack);
663  Int_t IdTrackMatching=-1;
664  //if(ctheta) (NADC*0.5*stheta + rho)/ctheta;
665  HarpoRecoHoughTracks* tr = new HarpoRecoHoughTracks(mtrack,theta,rho,sigTheta,sigRho,Qtrack,AngleTrack,x0,typeTrack,IdTrackMatching,Ncluster,plane,xStart,zStart,xEnd,zEnd);
666 
667  fEvt->GetRecoEvent()->AddTracks(tr);
668  tr->Delete();
669  if(gHarpoDebug>0)
670  Info("FindTrack","Found track %d",mtrack);
671  mtrack++;
672  } while(mtrack<NTRACK);
673 
674 }
675 
676 void HarpoHoughTracking::InitTrack()
677 {
678  int j;
679  // mtrack=0;
680  Qtrack=0;
681 
682  fNloopMax = -1;
683  fTrackIdMax = -1;
684  for(Int_t plane = 0; plane<2; plane++){
685  for(Int_t track = 0; track<5; track++){
686  for(Int_t loop = 0; loop<5; loop++){
687  fHist[plane][track][loop] = 0;
688  }
689  }
690  }
691  for(j=0;j<NCLUMAX;j++){
692  iclusttrack[j]=-1;
693  }
694 
695 
696  for(Int_t itr = 0; itr<NTRACK; itr++){
697  if(fGraphRoute[itr])
698  fGraphRoute[itr]->Delete();
699  fGraphRoute[itr] = new TGraph();
700  }
701 
702 
703 }
704 
705  void HarpoHoughTracking::Init()
706  {
707 
708  Info("Init","initializing tracking");
709 
710 
711  for(Int_t itr = 0; itr<NTRACK; itr++){
712  fGraphRoute[itr] = 0;
713  fGraphRoute[itr] = 0;
714  }
715 
716  // fCanvas = 0;
717  // fTracksCanvas = 0;
718  fDisplay = 0;
719  fSave = kFALSE;
720  fChooseQminC = 0;
721 
722 
723  // fSigmaMaxRebin = 0.25;
724  // fMinPairs = 50;
725  // fCutoffHough = 6.;
726  // RouteHoughTight = 3.;
727  // RouteHoughLoose = 6.;
728  // fDistMin = 50;
729  // fNclMin = 10;
730  fSigmaMaxRebin = 0.25;
731  fMinPairs = 25;
732  fCutoffHough = 6.;
733  RouteHoughTight = 3.;
734  RouteHoughLoose = 6.;
735  fDistMin = 1;
736  fNclMin = 5;
737 
738  fQmin[0] = 500;
739  fQmin[1] = 500;
740  fQmin[2] = 500;
741 
742  fWmin[0] = 10;
743  fWmin[1] = 3;
744  fWmin[2] = -1;
745  fWmax[0] = 30;
746  fWmax[1] = 30;
747  fWmax[2] = NADC*NCHAN;
748  // mtrack = 0;
749 
750  Long64_t WminT=80;
751  if( ! gHConfig->Lookup("hough.fWminT",WminT ))
752  Info("Constructor","Use default fWminT %d",fWmin[TCLUSTER]);
753  else
754  fWmin[TCLUSTER] = WminT;
755  Long64_t WminC=80;
756  if( ! gHConfig->Lookup("hough.fWminC",WminC ))
757  Info("Constructor","Use default fWminC %d",fWmin[CCLUSTER]);
758  else
759  fWmin[CCLUSTER] = WminC;
760 
761  Long64_t WmaxT=80;
762  if( ! gHConfig->Lookup("hough.fWmaxT",WmaxT ))
763  Info("Constructor","Use default fWmaxT %d",fWmax[TCLUSTER]);
764  else
765  fWmax[TCLUSTER] = WmaxT;
766  Long64_t WmaxC=80;
767  if( ! gHConfig->Lookup("hough.fWmaxC",WmaxC ))
768  Info("Constructor","Use default fWmaxC %d",fWmax[CCLUSTER]);
769  else
770  fWmax[CCLUSTER] = WmaxC;
771 
772  Long64_t QminT=80;
773  if( ! gHConfig->Lookup("hough.fQminT",QminT ))
774  Info("Constructor","Use default fQminT %d",fQmin[TCLUSTER]);
775  else
776  fQmin[TCLUSTER] = QminT;
777  Long64_t QminC=80;
778  if( ! gHConfig->Lookup("hough.fQminC",QminC ))
779  Info("Constructor","Use default fQminC %d",fQmin[CCLUSTER]);
780  else
781  fQmin[CCLUSTER] = QminC;
782  Long64_t QminB=80;
783  if( ! gHConfig->Lookup("hough.fQminBloc",QminB ))
784  Info("Constructor","Use default fQminBloc %d",fQmin[BLOCCLUSTER]);
785  else
786  fQmin[BLOCCLUSTER] = QminB;
787 
788  Long64_t MinPairs = 50;
789  if( ! gHConfig->Lookup("hough.fMinPairs",MinPairs ))
790  Info("Constructor","Use default fMinPairs %d",fMinPairs);
791  else
792  fMinPairs = MinPairs;
793 
794 
795  Double_t cutoffhough = 6., routehoughTight = 3., routehoughLoose = 6.;//, restcluster=20.;
796  if( ! gHConfig->Lookup("hough.CutoffHough", cutoffhough))
797  Info("Constructor","Use default CutoffHough %g",fCutoffHough);
798  else
799  fCutoffHough = cutoffhough;
800  if( ! gHConfig->Lookup("hough.RouteHoughTight", routehoughTight))
801  Info("Constructor","Use default RouteHoughTight %g",RouteHoughTight);
802  else
803  RouteHoughTight = routehoughTight;
804  if( ! gHConfig->Lookup("hough.RouteHoughLoose", routehoughLoose))
805  Info("Constructor","Use default RouteHoughLoose %g",RouteHoughLoose);
806  else
807  RouteHoughLoose = routehoughLoose;
808  Double_t sigmaMaxRebin = 0.25;
809  if( ! gHConfig->Lookup("hough.SigmaMaxRebin", sigmaMaxRebin))
810  Info("Constructor","Use default SigmaMaxRebin %g",fSigmaMaxRebin);
811  else
812  fSigmaMaxRebin = sigmaMaxRebin;
813  // if( ! gHConfig->Lookup("hough.RestCluster", restcluster))
814  // Info("Constructor","Use default RestCluster %g",RestCluster);
815  // else
816  // RestCluster = restcluster;
817 
818  // fNbinr = 32;
819  // fNbint = 80;
820  fNbinr = 128;
821  fNbint = 320;
822  Long64_t Nbinr=32, Nbint=80;
823  if( ! gHConfig->Lookup("hough.fNbinr",Nbinr ))
824  Info("Constructor","Use default fNbinr %d",fNbinr);
825  else
826  fNbinr = Nbinr;
827 
828  if( ! gHConfig->Lookup("hough.fNbint",Nbint ))
829  Info("Constructor","Use default fNbint %d",fNbint);
830  else
831  fNbint = Nbint;
832 
833 
834  hSigmaRho = new TH1F("hSigmaRho","",100,0,1);
835  hSigmaRho->SetLineColor(kRed);
836  hSigmaTheta = new TH1F("hSigmaTheta","",100,0,1);
837  hSigmaTheta->SetLineColor(kRed);
838  hSigmaRhoEvent = new TH1F("hSigmaRhoEvent","",100,0,1);
839  hSigmaThetaEvent = new TH1F("hSigmaThetaEvent","",100,0,1);
840  hSigmaTheta->SetLineColor(kRed);
841 
842  }
843 
844 void HarpoHoughTracking::Save(char * /* mode */)
845 {
846 
847  // OpenHistFile("recoHoughTracking");
848 
849 
850 }
851 
852 
853 
854 
855 
856 
857 
858 
859 
860 
861 
862 //TGFrame* HarpoHoughTracking::ConfigFrame(TGFrame* f, UInt_t xsize, UInt_t ysize, TGMainFrame* fMain, Int_t id)
863 void HarpoHoughTracking::ConfigFrame(TGMainFrame* fMain, Int_t id)
864 {
865 
866  UInt_t xsize = 200;
867  UInt_t ysize2 = 20;
868  UInt_t ysize = 9*ysize2;
869  TGTransientFrame* main = new TGTransientFrame(gClient->GetRoot(), fMain, xsize, ysize);
870  main->Connect("CloseWindow()", "HarpoRecoMonitorGui", main, "CloseWindow()");
871  main->DontCallClose(); // to avoid double deletions.
872 
873  // use hierarchical cleaning
874  main->SetCleanup(kDeepCleanup);
875 
876  TGVerticalFrame* f213 = new TGVerticalFrame(main,xsize,ysize,kVerticalFrame);
877  TGLabel* fTrackingLabel = new TGLabel(f213, "Hough Tracking");
878 
879  TGHorizontalFrame* wMinFrameC = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
880  TGLabel* wMinLabelC = new TGLabel(wMinFrameC, "WminCl");
881  fChooseWminC = new TGNumberEntry(wMinFrameC);
882  fChooseWminC->SetNumber(fWmin[CCLUSTER]);
883 
884  TGHorizontalFrame* wMaxFrameC = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
885  TGLabel* wMaxLabelC = new TGLabel(wMaxFrameC, "WmaxCl");
886  fChooseWmaxC = new TGNumberEntry(wMaxFrameC);
887  fChooseWmaxC->SetNumber(fWmax[CCLUSTER]);
888 
889  TGHorizontalFrame* qMinFrameC = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
890  TGLabel* qMinLabelC = new TGLabel(qMinFrameC, "QminCl");
891  fChooseQminC = new TGNumberEntry(qMinFrameC);
892  fChooseQminC->SetNumber(fQmin[CCLUSTER]);
893 
894 
895  TGHorizontalFrame* wMinFrameT = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
896  TGLabel* wMinLabelT = new TGLabel(wMinFrameT, "WminCl");
897  fChooseWminT = new TGNumberEntry(wMinFrameT);
898  fChooseWminT->SetNumber(fWmin[TCLUSTER]);
899 
900  TGHorizontalFrame* wMaxFrameT = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
901  TGLabel* wMaxLabelT = new TGLabel(wMaxFrameT, "WmaxCl");
902  fChooseWmaxT = new TGNumberEntry(wMaxFrameT);
903  fChooseWmaxT->SetNumber(fWmax[TCLUSTER]);
904 
905  TGHorizontalFrame* qMinFrameT = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
906  TGLabel* qMinLabelT = new TGLabel(qMinFrameT, "QminCl");
907  fChooseQminT = new TGNumberEntry(qMinFrameT);
908  fChooseQminT->SetNumber(fQmin[TCLUSTER]);
909 
910 
911 
912  TGHorizontalFrame* nMinPairsFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
913  TGLabel* nMinPairsLabel = new TGLabel(nMinPairsFrame, "min pairs");
914  fChooseNminPairs = new TGNumberEntry(nMinPairsFrame);
915  fChooseNminPairs->SetNumber(fMinPairs);
916 
917  TGHorizontalFrame* CutoffHoughFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
918  TGLabel* CutoffHoughLabel = new TGLabel(CutoffHoughFrame, "Cutoff");
919  fChooseCutoffHough = new TGNumberEntry(CutoffHoughFrame);
920  fChooseCutoffHough->SetNumber(fCutoffHough);
921 
922  TGHorizontalFrame* RouteHoughTightFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
923  TGLabel* RouteHoughTightLabel = new TGLabel(RouteHoughTightFrame, "Tight Route");
924  fChooseRouteHoughTight = new TGNumberEntry(RouteHoughTightFrame);
925  fChooseRouteHoughTight->SetNumber(RouteHoughTight);
926 
927  TGHorizontalFrame* RouteHoughLooseFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
928  TGLabel* RouteHoughLooseLabel = new TGLabel(RouteHoughLooseFrame, "Loose Route");
929  fChooseRouteHoughLoose = new TGNumberEntry(RouteHoughLooseFrame);
930  fChooseRouteHoughLoose->SetNumber(RouteHoughLoose);
931 
932  TGHorizontalFrame* SigmaMaxRebinFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
933  TGLabel* SigmaMaxRebinLabel = new TGLabel(SigmaMaxRebinFrame, "#sigma_{rebin}");
934  fChooseSigmaMaxRebin = new TGNumberEntry(SigmaMaxRebinFrame);
935  fChooseSigmaMaxRebin->SetNumber(fSigmaMaxRebin);
936 
937  TGHorizontalFrame* DistMinFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
938  TGLabel* DistMinLabel = new TGLabel(DistMinFrame, "min distance");
939  fChooseDistMin = new TGNumberEntry(DistMinFrame);
940  fChooseDistMin->SetNumber(fDistMin);
941 
942  TGHorizontalFrame* NclMinFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
943  TGLabel* NclMinLabel = new TGLabel(NclMinFrame, "min N_{cl}");
944  fChooseNclMin = new TGNumberEntry(NclMinFrame);
945  fChooseNclMin->SetNumber(fNclMin);
946 
947  TGHorizontalFrame* SaveFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
948  TGLabel* SaveLabel = new TGLabel(SaveFrame, "Save hist");
949  fCheckSave = new TGCheckButton(SaveFrame);
950  //fCheckSave->SetNumber(fNclMin);
951 
952 
953  TGTextButton* setConf = new TGTextButton(f213,"Save Config",id);
954  setConf->Associate(fMain);
955 
956 
957  TGLayoutHints* fLayout1 = new TGLayoutHints(kLHintsLeft ,5,5,5,5);
958  TGLayoutHints* fLayout2 = new TGLayoutHints(kLHintsRight ,5,5,5,5);
959  TGLayoutHints* fLayout3 = new TGLayoutHints(kLHintsTop | kLHintsExpandX ,5,5,5,5);
960  TGLayoutHints* fLayout4 = new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,5,5,5,5);
961 
962  main->AddFrame(f213,fLayout4);
963  f213->AddFrame(fTrackingLabel,fLayout3);
964 
965  f213->AddFrame(wMinFrameC,fLayout3);
966  wMinFrameC->AddFrame(wMinLabelC,fLayout1);
967  wMinFrameC->AddFrame(fChooseWminC,fLayout2);
968  f213->AddFrame(wMaxFrameC,fLayout3);
969  wMaxFrameC->AddFrame(wMaxLabelC,fLayout1);
970  wMaxFrameC->AddFrame(fChooseWmaxC,fLayout2);
971  f213->AddFrame(qMinFrameC,fLayout3);
972  qMinFrameC->AddFrame(qMinLabelC,fLayout1);
973  qMinFrameC->AddFrame(fChooseQminC,fLayout2);
974 
975  f213->AddFrame(wMinFrameT,fLayout3);
976  wMinFrameT->AddFrame(wMinLabelT,fLayout1);
977  wMinFrameT->AddFrame(fChooseWminT,fLayout2);
978  f213->AddFrame(wMaxFrameT,fLayout3);
979  wMaxFrameT->AddFrame(wMaxLabelT,fLayout1);
980  wMaxFrameT->AddFrame(fChooseWmaxT,fLayout2);
981  f213->AddFrame(qMinFrameT,fLayout3);
982  qMinFrameT->AddFrame(qMinLabelT,fLayout1);
983  qMinFrameT->AddFrame(fChooseQminT,fLayout2);
984 
985  f213->AddFrame(nMinPairsFrame,fLayout3);
986  nMinPairsFrame->AddFrame(nMinPairsLabel,fLayout1);
987  nMinPairsFrame->AddFrame(fChooseNminPairs,fLayout2);
988 
989  f213->AddFrame(CutoffHoughFrame,fLayout3);
990  CutoffHoughFrame->AddFrame(CutoffHoughLabel,fLayout1);
991  CutoffHoughFrame->AddFrame(fChooseCutoffHough,fLayout2);
992  f213->AddFrame(RouteHoughTightFrame,fLayout3);
993  RouteHoughTightFrame->AddFrame(RouteHoughTightLabel,fLayout1);
994  RouteHoughTightFrame->AddFrame(fChooseRouteHoughTight,fLayout2);
995  f213->AddFrame(RouteHoughLooseFrame,fLayout3);
996  RouteHoughLooseFrame->AddFrame(RouteHoughLooseLabel,fLayout1);
997  RouteHoughLooseFrame->AddFrame(fChooseRouteHoughLoose,fLayout2);
998  f213->AddFrame(DistMinFrame,fLayout3);
999  DistMinFrame->AddFrame(DistMinLabel,fLayout1);
1000  DistMinFrame->AddFrame(fChooseDistMin,fLayout2);
1001  f213->AddFrame(NclMinFrame,fLayout3);
1002  NclMinFrame->AddFrame(NclMinLabel,fLayout1);
1003  NclMinFrame->AddFrame(fChooseNclMin,fLayout2);
1004  f213->AddFrame(SigmaMaxRebinFrame,fLayout3);
1005  SigmaMaxRebinFrame->AddFrame(SigmaMaxRebinLabel,fLayout1);
1006  SigmaMaxRebinFrame->AddFrame(fChooseSigmaMaxRebin,fLayout2);
1007  f213->AddFrame(SaveFrame,fLayout3);
1008  SaveFrame->AddFrame(SaveLabel,fLayout1);
1009  SaveFrame->AddFrame(fCheckSave,fLayout2);
1010 
1011  f213->AddFrame(setConf,fLayout3);
1012 
1013 
1014  main->MapSubwindows();
1015  main->MapWindow();
1016  main->Resize();
1017  return;
1018 }
1019 
1020 
1021 void HarpoHoughTracking::DisplayAnalysis(TRootEmbeddedCanvas* ecTab, TGListBox* /* infobox */)
1022 {
1023 
1024  // Bool_t save = kTRUE;
1025  TFile* savefile = NULL;
1026  // if(fSave) savefile = new TFile("tracking.root","recreate");
1027 
1028  Char_t filename[256];
1029  sprintf(filename,"plots/tracking/run%d_evt%ld.root",fEvt->GetRunHeader()->GetRun(),fEvt->GetHeader()->GetEvtNo());
1030  if(fSave) savefile = new TFile(filename,"recreate");
1031 
1032  Info("DisplayAnalysis","%d %d",fNloopMax, fTrackIdMax);
1033  if(fNloopMax<0) return;
1034  if(fTrackIdMax<0) return;
1035 
1036  TCanvas* fCanvas = ecTab->GetCanvas();
1037  fCanvas->GetPad(1)->Delete();
1038  fCanvas->GetPad(2)->Delete();
1039  fCanvas->Divide(2);
1040  fCanvas->GetPad(1)->Divide(fNloopMax+1,fTrackIdMax+2);
1041  fCanvas->GetPad(2)->Divide(fNloopMax+1,fTrackIdMax+2);
1042  TLatex* latex = new TLatex();
1043  latex->SetTextFont(132);
1044  latex->SetTextAlign(12);
1045  latex->SetTextSize(0.08);
1046 
1047  for(Int_t plane = 0; plane<2; plane++){
1048  for(Int_t track = 0; track<=fTrackIdMax; track++){
1049  for(Int_t loop = 0; loop<=fNloopMax; loop++){
1050  Info("DisplayAnalysis","%d %d %d",plane,track,loop);
1051  if(fHist[plane][track][loop] != 0){
1052  if(fSave) fHist[plane][track][loop]->Write(Form("h_%d_%d_%d",plane,track,loop));
1053  MakeNiceHisto(fHist[plane][track][loop],fCanvas->GetPad(plane+1)->GetPad(track*(fNloopMax+1)+loop+1),"col");
1054  Double_t xmin = fHist[plane][track][loop]->GetXaxis()->GetXmin();
1055  Double_t xmax = fHist[plane][track][loop]->GetXaxis()->GetXmax();
1056  Double_t ymin = fHist[plane][track][loop]->GetYaxis()->GetXmin();
1057  if(ymin<-1.6) ymin = -TMath::Pi()/2;
1058  Double_t ymax = fHist[plane][track][loop]->GetYaxis()->GetXmax();
1059  if(ymax>1.6) ymax = TMath::Pi()/2;
1060  fHist[plane][track][loop]->Delete();
1061  latex->DrawLatex(xmin + 0.1*(xmax-xmin),ymin + 0.9*(ymax-ymin),Form("#sigma_{#rho} = %g",fSigmaRhoPeak[plane][track][loop]));
1062  latex->DrawLatex(xmin + 0.1*(xmax-xmin),ymin + 0.8*(ymax-ymin),Form("#sigma_{#theta} = %g",fSigmaThetaPeak[plane][track][loop]));
1063  latex->DrawLatex(xmin + 0.1*(xmax-xmin),ymin + 0.7*(ymax-ymin),Form("N_{pairs} = %g",fNpairsPeak[plane][track][loop]));
1064  }
1065  }
1066  }
1067  }
1068  if(fSave) savefile->Close();
1069  fCanvas->GetPad(1)->cd((fNloopMax+1)*(fTrackIdMax+2));
1070  MakeNiceHisto(hSigmaRho,fCanvas->GetPad(1)->GetPad((fNloopMax+1)*(fTrackIdMax+2)));
1071  hSigmaRhoEvent->Scale(hSigmaRho->GetEntries()*1./hSigmaRhoEvent->GetEntries());
1072  hSigmaRhoEvent->DrawCopy("same");
1073  MakeNiceHisto(hSigmaTheta,fCanvas->GetPad(2)->GetPad((fNloopMax+1)*(fTrackIdMax+2)));
1074  hSigmaThetaEvent->Scale(hSigmaTheta->GetEntries()*1./hSigmaThetaEvent->GetEntries());
1075  hSigmaThetaEvent->Draw("same");
1076  fCanvas->Update();
1077 }
1078 
1079 
1080 void HarpoHoughTracking::SetConfig()
1081 {
1082 
1083  Info("SetConfig","Renewing configuration");
1084 
1085  fDisplay = 1;
1086  Info("SetConfig","###################### (%d)",fDisplay);
1087 
1088  if(!fChooseQminC) return;
1089 
1090  fQmin[CCLUSTER] = fChooseQminC->GetNumber();
1091  fQmin[TCLUSTER] = fChooseQminT->GetNumber();
1092  fQmin[BLOCCLUSTER] = fChooseQminC->GetNumber();
1093  fWmin[CCLUSTER] = fChooseWminC->GetNumber();
1094  fWmax[TCLUSTER] = fChooseWmaxT->GetNumber();
1095  fWmin[CCLUSTER] = fChooseWminC->GetNumber();
1096  fWmax[TCLUSTER] = fChooseWmaxT->GetNumber();
1097  fMinPairs = fChooseNminPairs->GetNumber();
1098  fSave = fCheckSave->IsOn();
1099 
1100 
1101  fCutoffHough = fChooseCutoffHough->GetNumber();
1102  RouteHoughLoose = fChooseRouteHoughLoose->GetNumber();
1103  RouteHoughTight = fChooseRouteHoughTight->GetNumber();
1104  fSigmaMaxRebin = fChooseSigmaMaxRebin->GetNumber();
1105  fDistMin = fChooseDistMin->GetNumber();
1106  fNclMin = fChooseNclMin->GetNumber();
1107 
1108  // fNbinr = 32;
1109  // fNbint = 80;
1110 }
HarpoHoughTracking::kfWidthTheta
static const Double_t kfWidthTheta
Definition: HarpoHoughTracking.h:71
NTRACK
#define NTRACK
Definition: HarpoHoughTracking.h:22
HarpoHoughTracking::kfWidthRho
static const Double_t kfWidthRho
Definition: HarpoHoughTracking.h:70
NCHAN
#define NCHAN
Definition: HarpoDccMap.h:16
HarpoHoughTracking::FindPairNew
void FindPairNew(Int_t Plan, Int_t trackId, Double_t &rho, Double_t &theta, Double_t &sigRho, Double_t &sigTheta)
Ovreloaded medod whic do all job.
Definition: HarpoHoughTracking.cxx:109
HarpoRecoClusters::AddIdClusterTrack
void AddIdClusterTrack(Int_t val)
Definition: HarpoRecoEvent.cxx:75
main
int main(int argc, char **argv)
Definition: drawAnalyse.cxx:25
NCLUMAX
#define NCLUMAX
Definition: HarpoHoughTracking.h:23
HarpoDetSet::isExist
Bool_t isExist(ULong_t det)
Detecror date exist //! Number of Real Detectors.
Definition: HarpoDetSet.h:33
HarpoRecoEvent
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
Definition: HarpoRecoEvent.h:474
BLOCCLUSTER
#define BLOCCLUSTER
Definition: HarpoRecoEvent.h:24
HarpoDccMap::print
virtual void print()
Definition: HarpoDccMap.cxx:124
HarpoEventHeader
FullEvent Header not scecific to the detectors The class is ....
Definition: HarpoEventHeader.h:27
HarpoHoughTracking::DisplayAnalysis
void DisplayAnalysis(TRootEmbeddedCanvas *ecTab, TGListBox *info)
Definition: HarpoHoughTracking.cxx:1021
HarpoRecoClusters::RemoveAllClusterTrack
void RemoveAllClusterTrack()
Definition: HarpoRecoEvent.h:74
HarpoEventHeader::print
virtual void print()
Definition: HarpoEventHeader.cxx:21
HarpoRecoClusters
Cluster object, containing position, charge and quality information.
Definition: HarpoRecoEvent.h:28
HarpoHoughTracking::Save
void Save(char *mode=NULL)
Definition: HarpoHoughTracking.cxx:844
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
HarpoHoughTracking::FindTrack
void FindTrack(Int_t plane)
Definition: HarpoHoughTracking.cxx:506
HarpoHoughTracking::kfFacHough
static const Double_t kfFacHough
Definition: HarpoHoughTracking.h:72
HOUGH
#define HOUGH
Definition: HarpoRecoEvent.h:25
MakeNiceHisto
TVirtualPad * MakeNiceHisto(TH1 *hist, TVirtualPad *c1, const char *opt, Bool_t copy)
Definition: MakeNiceHisto.cxx:221
HarpoHoughTracking::ConfigFrame
void ConfigFrame(TGMainFrame *fMain, Int_t id)
Definition: HarpoHoughTracking.cxx:863
NADC
#define NADC
Definition: HarpoDccMap.h:18
TCLUSTER
#define TCLUSTER
Definition: HarpoRecoEvent.h:23
gHarpoDebug
Long64_t gHarpoDebug
Definition: HarpoDebug.cxx:9
HarpoConfigFile::Lookup
Bool_t Lookup(const char *path, Bool_t &val)
Lookup function for scalar values.
Definition: HarpoConfigFile.cxx:85
HarpoHoughTracking
A virtual class which define intrafece between HARPO Reader and Event Analysis code.
Definition: HarpoHoughTracking.h:30
gkNDetectors
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
CCLUSTER
#define CCLUSTER
Definition: HarpoRecoEvent.h:22
gHConfig
HarpoConfig * gHConfig
Definition: HarpoConfig.cxx:101
HarpoRecoHoughTracks
HarpoRecoTracks object, obtained with Hough tracking method.
Definition: HarpoRecoEvent.h:190
HarpoRecoEvent::GetClustersArray
TClonesArray * GetClustersArray()
Definition: HarpoRecoEvent.h:594
gHDetSet
R__EXTERN HarpoDetSet * gHDetSet
Definition: HarpoDetSet.h:71
Generated by   doxygen 1.8.6