31 #include "FairRootManager.h"
32 #include "FairRunAna.h"
33 #include "FairRuntimeDb.h"
35 #include "TClonesArray.h"
40 #include "TSpectrum2.h"
41 #include "TSpectrum.h"
42 #include "TStopwatch.h"
60 double value=(par[0] * x - y + par[1]) /
TMath::Sqrt(par[0] * par[0] + 1);
69 TMatrixT<Double_t> *objtofit = (TMatrixT<Double_t> *) gMinuit->GetObjectFit();
72 Int_t hitcounter = objtofit->GetNrows();
74 for (Int_t ihit = 0; ihit < hitcounter; ihit++)
76 double r_reco =
fit_distance(objtofit[0][ihit][0], objtofit[0][ihit][1], par);
77 double delta = (objtofit[0][ihit][2] - r_reco)/objtofit[0][ihit][3];
79 cout <<
"reco iso " << r_reco << endl;
80 cout <<
"drift " << objtofit[0][ihit][2] << endl;
81 cout <<
"delta " << delta << endl;
82 chi2 += delta * delta;
91 TMatrixT<Double_t> *objtofit = (TMatrixT<Double_t> *) gMinuit->GetObjectFit();
94 Int_t hitcounter = objtofit->GetNrows();
96 for (Int_t ihit = 0; ihit < hitcounter; ihit++)
99 double delta =
TMath::Sqrt((objtofit[0][ihit][0] - par[0])*(objtofit[0][ihit][0] - par[0])+(objtofit[0][ihit][1] - par[1])*(objtofit[0][ihit][1] - par[1])) - par[2]+ par[3] * objtofit[0][ihit][2];
100 if(objtofit[0][ihit][2] == 0) chi2 += (delta * delta * 12.);
101 else chi2 += (delta*delta)/(pow(objtofit[0][ihit][3],2));
115 minimizer.SetPrintLevel(-1);
119 TMatrixT<Double_t> fitvect;
122 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++)
125 if(hit == fRefHit)
continue;
129 fitvect.ResizeTo(nofhits, 4);
132 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++)
135 if(hit == fRefHit)
continue;
160 mrk->SetMarkerColor(kRed);
165 mrk2->SetMarkerColor(kRed);
173 if(nofhits != counter) fitvect.ResizeTo(counter, 4);
179 minimizer.DefineParameter(0,
"m", mstart, 0.1, -3000., 3000.);
180 minimizer.DefineParameter(1,
"q", qstart, 0.1, -3000., 3000.);
182 minimizer.SetObjectFit(&fitvect);
183 minimizer.SetPrintLevel();
185 minimizer.SetMaxIterations(500);
189 minimizer.GetParameter(0, results[0], errors[0]);
190 minimizer.GetParameter(1, results[1], errors[1]);
192 cout <<
"fitm: " << results[0] << endl;
193 cout <<
"fitq: " << results[1] << endl;
196 if(fitm == 0)
return kFALSE;
228 minimizer.SetPrintLevel(-1);
232 TMatrixT<Double_t> fitvect;
235 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++)
238 if(hit == fRefHit)
continue;
242 fitvect.ResizeTo(nofhits, 4);
245 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++)
248 if(hit == fRefHit)
continue;
271 mrk->SetMarkerColor(kRed);
278 if(nofhits != counter) fitvect.ResizeTo(counter, 4);
284 minimizer.DefineParameter(0,
"x", xstart, 0.1, -3000., 3000.);
285 minimizer.DefineParameter(1,
"y", ystart, 0.1, -3000., 3000.);
286 minimizer.DefineParameter(2,
"R", rstart, 0.1, 0., 3000.);
287 minimizer.DefineParameter(3,
"sign", 0, 1, -1, 1);
289 minimizer.SetObjectFit(&fitvect);
290 minimizer.SetPrintLevel();
292 minimizer.SetMaxIterations(500);
296 minimizer.GetParameter(0, results[0], errors[0]);
297 minimizer.GetParameter(1, results[1], errors[1]);
298 minimizer.GetParameter(2, results[2], errors[2]);
299 minimizer.GetParameter(3, results[3], errors[3]);
301 cout <<
"xc: " << results[0] << endl;
302 cout <<
"yc: " << results[1] << endl;
303 cout <<
"R: " << results[2] << endl;
304 cout <<
"sign: " << results[3] << endl;
334 PndTrkTrackFinder::PndTrkTrackFinder() : FairTask(
"secondary track finder", 0), fDisplayOn(kFALSE), fPersistence(kTRUE), fUseMVDPix(kTRUE), fUseMVDStr(kTRUE), fUseSTT(kTRUE), fUseSCIT(kTRUE), fUseGEM(kTRUE), fSecondary(kFALSE), fInitDone(kFALSE), fMvdPix_RealDistLimit(1000), fMvdStr_RealDistLimit(1000), fStt_RealDistLimit(1000), fMvdPix_ConfDistLimit(1000), fMvdStr_ConfDistLimit(1000), fStt_ConfDistLimit(1000), fUmin(-0.07), fUmax(0.07), fVmin(-0.07), fVmax(0.07), fRmin(-1.5), fRmax(1.5), fThetamin(0), fThetamax(180), fRefHit(NULL), fDelPrim(kFALSE), fNofPrimaries(0) {
343 PndTrkTrackFinder::PndTrkTrackFinder(
int verbose) : FairTask(
"secondary track finder", verbose), fDisplayOn(kFALSE), fPersistence(kTRUE), fUseMVDPix(kTRUE), fUseMVDStr(kTRUE), fUseSTT(kTRUE), fUseSCIT(kTRUE), fUseGEM(kTRUE), fSecondary(kFALSE), fInitDone(kFALSE), fMvdPix_RealDistLimit(1000), fMvdStr_RealDistLimit(1000), fStt_RealDistLimit(1000), fMvdPix_ConfDistLimit(1000), fMvdStr_ConfDistLimit(1000), fStt_ConfDistLimit(1000), fUmin(-0.07), fUmax(0.07), fVmin(-0.07), fVmax(0.07), fRmin(-1.5), fRmax(1.5), fThetamin(0), fThetamax(180), fRefHit(NULL), fDelPrim(kFALSE), fNofPrimaries(0) {
423 FairRootManager* ioman = FairRootManager::Instance();
425 cout <<
"-E- PndTrkTrackFinder::Init: "
426 <<
"RootManager not instantiated, return!" << endl;
435 cout <<
"-W- PndTrkTrackFinder::Init: "
436 <<
"No STTHit array, return!" << endl;
443 std::cout <<
"-W- PndTrkTrackFinder::Init: " <<
"No MVD Pixel hitArray, return!" << std::endl;
450 std::cout <<
"-W- PndTrkTrackFinder::Init: " <<
"No MVD Strip hitArray, return!" << std::endl;
457 std::cout <<
"-W- PndTrkTrackFinder::Init: " <<
"No SciT hitArray, return!" << std::endl;
464 std::cout <<
"-W- PndTrkTrackFinder::Init: " <<
"No GEM hitArray, return!" << std::endl;
484 display =
new TCanvas(
"display",
"display", 0, 0, 800, 800);
495 fTimer =
new TStopwatch();
508 fLineHisto =
new TH2F(
"fLineHisto",
"hl", 360, -360, 360, 500, -400, 400);
529 FairRuntimeDb*
rtdb = FairRunAna::Instance()->GetRuntimeDb();
550 std::map< int, std::vector< int > > det_to_hitids;
553 cout <<
fNofPrimaries <<
" --->" << det_to_hitids.size() << endl;
558 std::map< int, bool > primaries =
PrimaryCheck(FairRootManager::Instance()->GetBranchId(
fSttBranch), det_to_hitids);
559 for(
size_t ihit = 0; ihit < primaries.size(); ihit++) {
560 if(primaries[ihit] ==
true)
continue;
570 for(
size_t ihit = 0; ihit < primaries.size(); ihit++) {
571 if(primaries[ihit] ==
true)
continue;
581 for(
size_t ihit = 0; ihit < primaries.size(); ihit++) {
582 if(primaries[ihit] ==
true)
continue;
605 std::map< int, bool > primaries =
PrimaryCheck(FairRootManager::Instance()->GetBranchId(
fGemBranch), det_to_hitids);
606 for(
size_t ihit = 0; ihit < primaries.size(); ihit++) {
607 if(primaries[ihit] ==
true)
continue;
634 if(
fVerbose > 0) cout <<
"PndTrkTrackFinder:: *********************** " <<
fEventCounter <<
" ***********************" << endl;
641 cout <<
"number of stt hits " <<
fSttHitArray->GetEntriesFast() << endl;
644 cout <<
"number of scit hits " <<
fSciTHitArray->GetEntriesFast() << endl;
645 cout <<
"number of gem hits " <<
fGemHitArray->GetEntriesFast() << endl;
663 cout <<
" STARTING" << endl;
689 double numx, numy, den;
699 int nind = indiv.GetEntriesFast();
700 if(nind == 0)
continue;
702 TArrayI indhitids(nind + 1);
705 for(
int jhit = 0; jhit < nind; jhit++) {
719 TVector3 indpos(numx, numy, 0.);
731 indhit->
Draw(kMagenta);
743 std::map< int, std::vector< int > > maplay2hits;
749 std::map< int, std::vector< int > >::iterator it = maplay2hits.begin();
750 it = maplay2hits.find(layerID);
751 if(it == maplay2hits.end()) {
752 std::vector< int >
hits;
753 hits.push_back(ihit);
754 maplay2hits[layerID] =
hits;
756 std::vector< int > hits2 = maplay2hits[layerID];
759 std::vector< int >
hits = maplay2hits[layerID];
760 hits.push_back(ihit);
761 maplay2hits[layerID] =
hits;
766 for(
size_t ilay = 0; ilay < maplay2hits.size(); ilay++) {
767 std::vector< int >
hits = maplay2hits[ilay];
771 std::vector< std::vector < int > > trackcandidates;
778 std::vector< int > hits0 = maplay2hits[0];
779 for(
size_t ihit = 0; ihit < hits0.size(); ihit++) {
780 std::vector< int > hits7 = maplay2hits[7];
781 for(
size_t jhit = 0; jhit < hits7.size(); jhit++) {
795 double cosalpha = hitipos * hitjpos / (hitipos.Mod() * hitjpos.Mod());
796 if(cosalpha < 0.94) {
801 std::vector< int > couple;
802 couple.push_back(hits0[ihit]);
803 couple.push_back(hits7[jhit]);
804 trackcandidates.push_back(couple);
829 std::vector< std::vector < int > > trackcandidates2;
830 for(
size_t ipair = 0; ipair < trackcandidates.size(); ipair++) {
831 std::vector< int > couple = trackcandidates[ipair];
833 std::vector< int > hits16 = maplay2hits[16];
834 for(
size_t jhit = 0; jhit < hits16.size(); jhit++) {
849 double cosalpha = hit1_0pos * hit1_16pos/ (hit1_0pos.Mod() * hit1_16pos.Mod());
851 if(cosalpha < 0.900) {
858 std::vector< int > triplet;
859 triplet.push_back(couple[0]);
860 triplet.push_back(couple[1]);
861 triplet.push_back(hits16[jhit]);
863 trackcandidates2.push_back(triplet);
889 std::vector< std::vector< int > > trackcandidates3;
890 for(
size_t iqua = 0; iqua < trackcandidates2.size(); iqua++) {
891 std::vector< int > triplet = trackcandidates2[iqua];
898 std::vector< int > hits20 = maplay2hits[20];
899 for(
size_t jhit = 0; jhit < hits20.size(); jhit++) {
911 double cosalpha = hit2_1pos * hit2_20pos/ (hit2_1pos.Mod() * hit2_20pos.Mod());
912 if(cosalpha < 0.94) {
919 std::vector< int > quadriplet;
920 quadriplet.push_back(triplet[0]);
921 quadriplet.push_back(triplet[1]);
922 quadriplet.push_back(triplet[2]);
923 quadriplet.push_back(hits20[jhit]);
925 trackcandidates3.push_back(quadriplet);
957 std::vector < std::vector < double > > tracks;
958 std::vector< std::vector< int > > trackcandidates4;
959 for(
size_t iqua = 0; iqua < trackcandidates3.size(); iqua++) {
960 std::vector< int > quadriplet = trackcandidates3[iqua];
966 double x01, y01, rad1;
968 double x02, y02, rad2;
970 double x03, y03, rad3;
973 double xm = 0, ym = 0, rm = 0;
974 int positive[2] = {0, 0}, negative[2] = {0, 0};
975 if(x01 > 0) positive[0]++;
977 if(x02 > 0) positive[0]++;
979 if(x03 > 0) positive[0]++;
982 if(y01 > 0) positive[1]++;
984 if(y02 > 0) positive[1]++;
986 if(y03 > 0) positive[1]++;
989 bool posit[2] = {
false,
false};
990 positive[0] > negative[0] ? posit[0] =
true : posit[0] =
false;
991 positive[1] > negative[1] ? posit[1] =
true : posit[1] =
false;
994 if(posit[0] ==
true && x01 > 0) {
995 if(posit[1] ==
true && y01 > 0) { xm += x01; ym += y01; rm += rad1; count++; }
996 else if(posit[1] ==
false && y01 < 0) { xm += x01; ym += y01; rm += rad1; count++; }
998 else if(posit[0] ==
false && x01 < 0) {
999 if(posit[1] ==
true && y01 > 0) { xm += x01; ym += y01; rm += rad1; count++; }
1000 else if(posit[1] ==
false && y01 < 0) { xm += x01; ym += y01; rm += rad1; count++; }
1003 if(posit[0] ==
true && x02 > 0) {
1004 if(posit[1] ==
true && y02 > 0) { xm += x02; ym += y02; rm += rad2; count++; }
1005 else if(posit[1] ==
false && y02 < 0) { xm += x02; ym += y02; rm += rad2; count++; }
1007 else if(posit[0] ==
false && x02 < 0) {
1008 if(posit[1] ==
true && y02 > 0) { xm += x02; ym += y02; rm += rad2; count++; }
1009 else if(posit[1] ==
false && y02 < 0) { xm += x02; ym += y02; rm += rad2; count++; }
1012 if(posit[0] ==
true && x03 > 0) {
1013 if(posit[1] ==
true && y03 > 0) { xm += x03; ym += y03; rm += rad3; count++; }
1014 else if(posit[1] ==
false && y03 < 0) { xm += x03; ym += y03; rm += rad3; count++; }
1016 else if(posit[0] ==
false && x03 < 0) {
1017 if(posit[1] ==
true && y03 > 0) { xm += x03; ym += y03; rm += rad3; count++; }
1018 else if(posit[1] ==
false && y03 < 0) { xm += x03; ym += y03; rm += rad3; count++; }
1033 double maxlimit = 1.;
1041 cout <<
"circ1 " << x01 <<
" " << y01 <<
" " << rad1 << endl;
1042 cout <<
"circ2 " << x02 <<
" " << y02 <<
" " << rad2 << endl;
1043 cout <<
"circ3 " << x03 <<
" " << y03 <<
" " << rad3 << endl;
1044 cout <<
"circm " << xm <<
" " << ym <<
" " << rm << endl;
1046 TArc *arc01 =
new TArc(x01, y01, rad1);
1047 arc01->SetFillStyle(0);
1048 arc01->SetLineColor(1);
1049 arc01->Draw(
"SAME");
1050 TArc *arc02 =
new TArc(x02, y02, rad2);
1051 arc02->SetFillStyle(0);
1052 arc02->SetLineColor(2);
1053 arc02->Draw(
"SAME");
1054 TArc *arc03 =
new TArc(x03, y03, rad3);
1055 arc03->SetFillStyle(0);
1056 arc03->SetLineColor(3);
1057 arc03->Draw(
"SAME");
1060 mrk1->SetMarkerColor(kBlack);
1064 mrk2->SetMarkerColor(kRed);
1068 mrk3->SetMarkerColor(kGreen);
1072 mrk4->SetMarkerColor(kBlue);
1075 TArc *arcm =
new TArc(xm, ym, rm);
1076 arcm->SetFillStyle(0);
1077 arcm->SetLineStyle(2);
1078 arcm->SetLineColor(kMagenta);
1088 cout <<
"next quadriplet?" << endl;
1093 std::vector< double >
track;
1094 track.push_back(xm);
1095 track.push_back(ym);
1096 track.push_back(rm);
1099 std::vector< std::vector < double > >::iterator itr = tracks.begin();
1100 itr = std::find(tracks.begin(), tracks.end(),
track);
1101 if(itr == tracks.end()) {
1108 if(rm > 10 && rm < 2500) {
1109 tracks.push_back(track);
1110 trackcandidates4.push_back(quadriplet);
1120 for(
size_t itrk = 0; itrk < tracks.size(); itrk++) {
1121 std::vector< double >
track = tracks[itrk];
1122 double x = track[0];
1123 double y = track[1];
1124 double r = track[2];
1141 std::vector < int > quadriplet = trackcandidates4[itrk];
1151 std::map< int, int > sectorids;
1158 if(sectorids.find(sec[
isec]) == sectorids.end()) sectorids[sec[isec]] = 1;
1159 else sectorids[sec[
isec]]++;
1164 int tmpsecentries = 0, tmpsec = -1;
1165 for(
size_t isec = 0;
isec < sectorids.size();
isec++) {
1166 if(tmpsecentries < sectorids[
isec]) {
1167 tmpsecentries = sectorids[
isec];
1171 int sectorID = tmpsec;
1175 bool border =
false;
1176 int othersecID = -1;
1177 for(
size_t isec = 0;
isec < sectorids.size();
isec++) {
1178 if(sectorids[
isec] > 0 && (
int)
isec != sectorID) {
1200 if(hit == hit0 || hit == hit1 || hit == hit2 || hit == hit3)
continue;
1203 if(border ==
false && tube->
GetSectorID() != sectorID)
continue;
1204 else if(border ==
true && (tube->
GetSectorID() != sectorID && tube->
GetSectorID() != othersecID))
continue;
1206 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(x, y)).Mod();
1207 double recoiso =
fabs(distance_hit_center - r);
1209 if(recoiso < 1.) cluster.
AddHit(hit);
1214 cluster.
Draw(kBlue);
1218 cout <<
"want to go on?" << endl;
1232 if(border ==
false && indhit->
GetSector() != sectorID)
continue;
1233 else if(border ==
true && (indhit->
GetSector() != sectorID && indhit->
GetSector() != othersecID))
continue;
1235 double distance_hit_center = (indhit->
GetPosition().XYvector() - TVector2(x, y)).Mod();
1236 double recoiso =
fabs(distance_hit_center - r);
1238 indcluster.
AddHit(indhit);
1255 double tmpdistance = -1;
1257 for(
int ihit = 0; ihit < indcluster.
GetNofHits(); ihit++) {
1260 if(position.Perp() > tmpdistance) {
1261 tmpdistance = position.Perp();
1266 if(tmphitID == -1)
continue;
1270 Double_t delta = 0, trasl[2] = {0., 0.};
1283 double rc_of_min, rc_of_max;
1287 double u = chit->
GetU();
1288 double v = chit->
GetV();
1303 double rimin, rimax;
1304 r1 < r2 ? (rimin =
r1, rimax =
r2) : (rimin = r2, rimax = r1);
1306 rimin <
fRmin ? (rc_of_min = rc,
fRmin = rimin) : fRmin;
1307 rimax >
fRmax ? (rc_of_max = rc,
fRmax = rimax) : fRmax;
1316 double delt =
fabs(dv - du)/2.;
1317 du < dv ? (fUmin -= delt,
fUmax += delt) : (fVmin -= delt,
fVmax += delt);
1322 for(
int ihit = 0; ihit < nofconfhits; ihit++) {
1343 std::vector< double > theta_max, r_max;
1344 std::vector < int > content_max;
1348 double fitm2, fitq2;
1352 TLine *line =
new TLine(-10.07, fitq2 + fitm2 * (-10.07), 10.07, fitq2 + fitm2 * (10.07));
1366 bool goodtrack =
true;
1367 double distance =
TMath::Sqrt((x - xc) * (x - xc) + (y - yc) * (y - yc));
1368 double r_minus_rc =
fabs(r - R);
1372 if(
fabs(distance - r_minus_rc) > 2.) {
1373 double recoisosum = 0;
1374 for(
int ihit = 0; ihit < cluster.
GetNofHits(); ihit++) {
1375 hit = cluster.
GetHit(ihit);
1376 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
1377 double recoiso =
fabs(distance_hit_center - R);
1378 recoisosum += recoiso;
1381 if(recoisosum > 10) goodtrack =
false;
1386 if(goodtrack ==
false)
continue;
1391 TArc *arcm =
new TArc(xc, yc, R);
1392 arcm->SetFillStyle(0);
1393 arcm->SetLineColor(kBlue);
1411 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
1412 double recoiso =
fabs(distance_hit_center - R);
1433 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
1434 double recoiso =
fabs(distance_hit_center - R);
1454 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
1455 double recoiso =
fabs(distance_hit_center - R);
1484 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
1485 double recoiso =
fabs(distance_hit_center - R);
1488 if(recoiso < tmpdistance) {
1490 tmpdistance = recoiso;
1500 cout <<
"herh" << endl;
1575 for(
int ihit = 0; ihit < nofconfhits; ihit++) {
1585 double xc2, yc2, R2;
1589 bool fitting =
AnalyticalFit(&cluster, xc, yc, R, fitm, fitq);
1590 if(fitting == kFALSE)
continue;
1594 double fita, fitb, fitc, epsilon;
1596 if(fitting == kFALSE)
continue;
1600 bool fitting =
MinuitFit(&cluster, fitm2, fitq2, fitm, fitq);
1601 if(fitting == kFALSE)
continue;
1608 TArc *arcm =
new TArc(xc2, yc2, R2);
1609 arcm->SetFillStyle(0);
1610 arcm->SetLineColor(5);
1618 distance =
TMath::Sqrt((xc - xc2) * (xc - xc2) + (yc - yc2) * (yc - yc2));
1619 r_minus_rc =
fabs(R - R2);
1623 if(
fabs(distance - r_minus_rc) > 2.)
continue;
1627 std::vector< std::pair< int, int > > dontuse;
1638 if(border ==
false && tube->
GetSectorID() != sectorID)
continue;
1639 else if(border ==
true && (tube->
GetSectorID() != sectorID && tube->
GetSectorID() != othersecID))
continue;
1640 double distance_hit_center = (tube->
GetPosition().XYvector() - TVector2(xc2, yc2)).Mod();
1641 double recoiso =
fabs(distance_hit_center - R2);
1653 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc2, yc2)).Mod();
1654 double recoiso =
fabs(distance_hit_center - R2);
1657 int samesensor =
false;
1663 double distance_hitk_center = (hitk->
GetPosition().XYvector() - TVector2(xc2, yc2)).Mod();
1664 double recoisok =
fabs(distance_hitk_center - R2);
1665 if(recoiso < recoisok) {
1668 dontuse.push_back(dontpair);
1673 dontuse.push_back(dontpair);
1687 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc2, yc2)).Mod();
1688 double recoiso =
fabs(distance_hit_center - R2);
1691 int samesensor =
false;
1697 double distance_hitk_center = (hitk->
GetPosition().XYvector() - TVector2(xc2, yc2)).Mod();
1698 double recoisok =
fabs(distance_hitk_center - R2);
1699 if(recoiso < recoisok) {
1702 dontuse.push_back(dontpair);
1707 dontuse.push_back(dontpair);
1720 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc2, yc2)).Mod();
1721 double recoiso =
fabs(distance_hit_center - R2);
1734 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc2, yc2)).Mod();
1735 double recoiso =
fabs(distance_hit_center - R2);
1737 if(recoiso < tmpdistance) {
1739 tmpdistance = recoiso;
1784 finaltrack.
Draw(kGreen);
1806 if(hit->
IsSttSkew() == kFALSE)
continue;
1808 if(border ==
false && hit->
GetSector() != sectorID)
continue;
1809 else if(border ==
true && (hit->
GetSector() != sectorID && hit->
GetSector() != othersecID))
continue;
1817 TVector3 first = tube->
GetPosition() + wireDirection * halflength;
1818 TVector3 second = tube->
GetPosition() - wireDirection * halflength;
1824 TVector2 intersection1, intersection2;
1827 if(nofintersections == 0)
continue;
1828 if(nofintersections >= 2) {
1829 if(
fVerbose) cout <<
"ERROR: MORE THAN 1 INTERSECTION!!" << endl;
1840 double beta = wireDirection.Phi();
1845 TVector2 rottangent(rtx, rty);
1846 rottangent = rottangent.Unit();
1857 Double_t rotm = rottangent.Y()/rottangent.X();
1867 x0a = (-rotp +
TMath::Sqrt(b * b + a * a * rotm * rotm)) / rotm;
1868 x0b = (-rotp -
TMath::Sqrt(b * b + a * a * rotm * rotm)) / rotm;
1871 double intxa = (x0a * b * b - rotm * rotp * a *
a) / (b * b + rotm * rotm * a * a);
1872 double intya = rotm * intxa + rotp;
1873 double intxb = (x0b * b * b - rotm * rotp * a *
a) / (b * b + rotm * rotm * a * a);
1874 double intyb = rotm * intxb + rotp;
1901 double y0a = y0anew;
1903 double y0b = y0bnew;
1909 TEllipse *ell1 =
new TEllipse(x0a, y0a, a, b, 0, 360, -beta);
1910 ell1->SetFillStyle(0);
1911 ell1->SetLineColor(4);
1913 TEllipse *ell2 =
new TEllipse(x0b, y0b, a, b, 0, 360, -beta);
1914 ell2->SetFillStyle(0);
1915 ell2->SetLineColor(6);
1918 TMarker *mrkinta =
new TMarker(intxa, intya, 20);
1919 mrkinta->SetMarkerColor(4);
1920 mrkinta->Draw(
"SAME");
1921 TMarker *mrkintb =
new TMarker(intxb, intyb, 20);
1922 mrkintb->SetMarkerColor(6);
1923 mrkintb->Draw(
"SAME");
1929 Double_t t = ((x0a + y0a) - (first.X() + first.Y())) / ((second.X() - first.X()) + (second.Y() - first.Y()));
1930 Double_t z0a = first.Z() + (second.Z() - first.Z()) * t;
1933 t = ((x0b + y0b) - (first.X() + first.Y())) / ((second.X() - first.X()) + (second.Y() - first.Y()));
1934 Double_t z0b = first.Z() + (second.Z() - first.Z()) * t;
1936 TVector3 center1(x0a, y0a, z0a);
1937 TVector3 center2(x0b, y0b, z0b);
1940 double dx = intxa - x0a;
1941 double dy = intya - y0a;
1942 TVector3 dxdy(dx, dy, 0.0);
1944 TVector3 tfirst = first + dxdy;
1945 TVector3 tsecond = second + dxdy;
1947 t = ((intxa + intya) - (tfirst.X() + tfirst.Y())) / ((tsecond.X() - tfirst.X()) + (tsecond.Y() - tfirst.Y()));
1948 double intza = tfirst.Z() + (tsecond.Z() - tfirst.Z()) * t;
1963 tfirst = first - dxdy;
1964 tsecond = second - dxdy;
1966 t = ((intxb + intyb) - (tfirst.X() + tfirst.Y())) / ((tsecond.X() - tfirst.X()) + (tsecond.Y() - tfirst.Y()));
1967 double intzb = tfirst.Z() + (tsecond.Z() - tfirst.Z()) * t;
1969 TVector3 fin_intersection1(intxa, intya, intza);
1970 TVector3 fin_intersection2(intxb, intyb, intzb);
1974 double phi1 = finaltrack.
ComputePhi(fin_intersection1);
1975 double phi2 = finaltrack.
ComputePhi(fin_intersection2);
1982 skewhitlist.
AddHit(
PndTrkSkewHit(*hit, trackID, center1, fin_intersection1, phi1, center2, fin_intersection2, phi2, a, b, -1, beta));
2087 int phisec[2] = {0, 0};
2088 for(
int ihit = 0; ihit < skewhitlist.
GetNofHits(); ihit++)
2093 TVector3 fin_intersection21(-999, -999, -999), fin_intersection22(-999, -999, -999);
2094 double phi21 = -999, phi22 = -999;
2096 if(hitj->
IsStt() == kTRUE) {
2104 double phimean = 0.5 * (phi21 + phi22);
2105 if(phimean >= 0 && phimean < 180) phisec[0]++;
2112 if(phi21 >= 0 && phi21 < 180) phisec[0]++;
2120 for(
int ihit = 0; ihit < skewhitlist.
GetNofHits(); ihit++)
2125 TVector3 fin_intersection21(-999, -999, -999), fin_intersection22(-999, -999, -999);
2126 double phi21 = -999, phi22 = -999;
2128 if(hitj->
IsStt() == kTRUE) {
2136 double phimean = 0.5 * (phi21 + phi22);
2138 if((phisec[0] > phisec[1]) && (phimean >= 180)) {
2144 else if((phisec[0] < phisec[1]) && (phimean < 180)) {
2156 if((phisec[0] > phisec[1]) && (phi21 >= 180)) {
2160 else if((phisec[0] < phisec[1]) && (phi21 < 180)) {
2170 std::vector < int > first, second;
2171 double fitm3 = 0, fitq3 = 0;
2172 for(
int ihit = 0; ihit < skewhitlist.
GetNofHits(); ihit++) {
2173 hit = skewhitlist.
GetHit(ihit);
2175 if(hit->
IsStt() == kFALSE)
continue;
2181 double phi1 = skewhit->
GetPhi1();
2182 double phi2 = skewhit->
GetPhi2();
2187 TLine *linezphi =
new TLine(phi1, fin_intersection1.Z(), phi2, fin_intersection2.Z());
2188 linezphi->SetLineStyle(1);
2189 linezphi->Draw(
"SAME");
2190 TMarker *mrkzphi1 =
new TMarker(phi1, fin_intersection1.Z(), 20);
2191 mrkzphi1->SetMarkerColor(kBlue - 9);
2192 mrkzphi1->Draw(
"SAME");
2193 TMarker *mrkzphi2 =
new TMarker(phi2, fin_intersection2.Z(), 20);
2194 mrkzphi2->SetMarkerColor(kMagenta - 7);
2195 mrkzphi2->Draw(
"SAME");
2269 for(
int ihit = 0; ihit < skewhitlist.
GetNofHits(); ihit++)
2271 hit = skewhitlist.
GetHit(ihit);
2274 TVector3 fin_intersection11(-999, -999, -999), fin_intersection12(-999, -999, -999);
2275 double phi11 = -999, phi12 = -999;
2277 if(hit->
IsStt() == kTRUE) {
2293 TLine *linezphi =
new TLine(phi11, fin_intersection11.Z(), phi12, fin_intersection12.Z());
2294 linezphi->SetLineStyle(1);
2295 linezphi->Draw(
"SAME");
2296 TMarker *mrkzphi1 =
new TMarker(phi11, fin_intersection11.Z(), 20);
2297 mrkzphi1->SetMarkerColor(kBlue - 9);
2298 mrkzphi1->Draw(
"SAME");
2299 TMarker *mrkzphi2 =
new TMarker(phi12, fin_intersection12.Z(), 20);
2300 mrkzphi2->SetMarkerColor(kMagenta - 7);
2301 mrkzphi2->Draw(
"SAME");
2304 TMarker *mrkzphi1 =
new TMarker(phi11, fin_intersection11.Z(), 20);
2305 mrkzphi1->SetMarkerColor(kGreen);
2306 mrkzphi1->Draw(
"SAME");
2312 for(
int jhit = ihit + 1; jhit < skewhitlist.
GetNofHits(); jhit++)
2317 TVector3 fin_intersection21(-999, -999, -999), fin_intersection22(-999, -999, -999);
2318 double phi21 = -999, phi22 = -999;
2320 if(hitj->
IsStt() == kTRUE) {
2335 double cost = (fin_intersection21.Z() - fin_intersection11.Z()) /
TMath::Sqrt((phi11 - phi21) * (phi11 - phi21) + (fin_intersection21.Z() - fin_intersection11.Z()) * (fin_intersection21.Z() - fin_intersection11.Z()));
2336 double theta = TMath::ACos(cost);
2337 double r1 = phi11 * cost + fin_intersection11.Z() *
TMath::Sin(theta);
2338 double r2 = phi21 * cost + fin_intersection21.Z() *
TMath::Sin(theta);
2339 if(
fabs(r1 - r2) > 1.e-9) {
2340 theta = -TMath::ACos(cost);
2341 r1 = phi11 * cost + fin_intersection11.Z() *
TMath::Sin(theta);
2343 fLineHisto->Fill(theta * TMath::RadToDeg(), r1);
2352 cost = (fin_intersection21.Z() - fin_intersection12.Z()) /
TMath::Sqrt((phi12 - phi21) * (phi12 - phi21) + (fin_intersection21.Z() - fin_intersection12.Z()) * (fin_intersection21.Z() - fin_intersection12.Z()));
2353 theta = TMath::ACos(cost);
2354 r1 = phi12 * cost + fin_intersection12.Z() *
TMath::Sin(theta);
2355 r2 = phi21 * cost + fin_intersection21.Z() *
TMath::Sin(theta);
2356 if(
fabs(r1 - r2) > 1.e-9) {
2357 theta = -TMath::ACos(cost);
2358 r1 = phi12 * cost + fin_intersection12.Z() *
TMath::Sin(theta);
2360 fLineHisto->Fill(theta * TMath::RadToDeg(), r1);
2370 cost = (fin_intersection11.Z() - fin_intersection22.Z()) /
TMath::Sqrt((phi22 - phi11) * (phi22 - phi11) + (fin_intersection11.Z() - fin_intersection22.Z()) * (fin_intersection11.Z() - fin_intersection22.Z()));
2371 theta = TMath::ACos(cost);
2372 r1 = phi22 * cost + fin_intersection22.Z() *
TMath::Sin(theta);
2373 r2 = phi11 * cost + fin_intersection11.Z() *
TMath::Sin(theta);
2374 if(
fabs(r1 - r2) > 1.e-9) {
2375 theta = -TMath::ACos(cost);
2376 r1 = phi22 * cost + fin_intersection22.Z() *
TMath::Sin(theta);
2378 fLineHisto->Fill(theta * TMath::RadToDeg(), r1);
2386 if(phi12 != -999 && phi22 != -999) {
2388 cost = (fin_intersection12.Z() - fin_intersection22.Z()) /
TMath::Sqrt((phi22 - phi12) * (phi22 - phi12) + (fin_intersection12.Z() - fin_intersection22.Z()) * (fin_intersection12.Z() - fin_intersection22.Z()));
2389 theta = TMath::ACos(cost);
2390 r1 = phi22 * cost + fin_intersection22.Z() *
TMath::Sin(theta);
2391 r2 = phi12 * cost + fin_intersection12.Z() *
TMath::Sin(theta);
2392 if(
fabs(r1 - r2) > 1.e-9) {
2393 theta = -TMath::ACos(cost);
2394 r1 = phi22 * cost + fin_intersection22.Z() *
TMath::Sin(theta);
2396 fLineHisto->Fill(theta * TMath::RadToDeg(), r1);
2417 int binx, biny, binz;
2418 fLineHisto->GetBinXYZ(bin, binx, biny, binz);
2419 double tpeak =
fLineHisto->GetXaxis()->GetBinCenter(binx);
2420 double rpeak =
fLineHisto->GetYaxis()->GetBinCenter(biny);
2423 fitq3 = rpeak/
TMath::Sin(tpeak * TMath::DegToRad());
2428 for(
int ihit = 0; ihit < skewhitlist.
GetNofHits(); ihit++)
2430 hit = skewhitlist.
GetHit(ihit);
2433 TVector3 fin_intersection11(-999, -999, -999), fin_intersection12(-999, -999, -999);
2434 double phi11 = -999, phi12 = -999;
2436 double phi = -999,
z = -999,
sigma = -999;
2438 if(hit->
IsStt() == kTRUE) {
2446 z = 0.5 * (fin_intersection11.Z() + fin_intersection12.Z());
2447 phi = 0.5 * (phi11 + phi12);
2448 sigma = (fin_intersection11.Z() - fin_intersection12.Z())/
TMath::Sqrt(12.);
2454 z = fin_intersection11.Z();
2463 TLine *linezphi =
new TLine(phi11, fin_intersection11.Z(), phi12, fin_intersection12.Z());
2464 linezphi->SetLineStyle(1);
2465 linezphi->Draw(
"SAME");
2466 TMarker *mrkzphi1 =
new TMarker(phi11, fin_intersection11.Z(), 20);
2467 mrkzphi1->SetMarkerColor(kBlue - 9);
2468 mrkzphi1->Draw(
"SAME");
2469 TMarker *mrkzphi2 =
new TMarker(phi12, fin_intersection12.Z(), 20);
2470 mrkzphi2->SetMarkerColor(kMagenta - 7);
2471 mrkzphi2->Draw(
"SAME");
2474 TMarker *mrkzphi1 =
new TMarker(phi11, fin_intersection11.Z(), 20);
2475 mrkzphi1->SetMarkerColor(kGreen);
2476 mrkzphi1->Draw(
"SAME");
2479 TMarker *mrkzphi =
new TMarker(phi,
z, 20);
2480 mrkzphi->SetMarkerColor(kOrange);
2481 mrkzphi->Draw(
"SAME");
2487 double fitm3b, fitq3b;
2493 if(
fabs(fitm3b) > 1e-10) {
2494 if(
fabs(fitm3b) > 0.1) {
2502 TLine *line =
new TLine(0, fitq3, 360, 360 * fitm3 + fitq3);
2503 line->SetLineColor(kOrange);
2514 for(
int ihit = 0; ihit < skewhitlist.
GetNofHits(); ihit++) {
2521 std::vector < std::pair< int, int > >::iterator itr = dontuse.begin();
2522 itr = std::find(dontuse.begin(), dontuse.end(), thispair);
2523 if(itr != dontuse.end())
continue;
2526 if(zdistance < 5) skewhitlist2.
AddHit(hit);
2527 else dontuse.push_back(thispair);
2530 TVector3 fin_intersection1 = ((
PndTrkSkewHit*) hit)->GetIntersection1();
2531 TVector3 fin_intersection2 = ((
PndTrkSkewHit*) hit)->GetIntersection2();
2535 double dist1 =
fabs(fitm3 * phi1 - fin_intersection1.Z() + fitq3) /
TMath::Sqrt(fitm3 * fitm3 + 1);
2536 double dist2 =
fabs(fitm3 * phi2 - fin_intersection2.Z() + fitq3) /
TMath::Sqrt(fitm3 * fitm3 + 1);
2539 dist1 < dist2 ? (distance = dist1, hit->
SetPosition(fin_intersection1), hit->
SetPhi(phi1)) : (distance = dist2, hit->
SetPosition(fin_intersection2), hit->
SetPhi(phi2));
2541 if(distance < 10) skewhitlist2.
AddHit(hit);
2547 double mvdmeandist = 0, sttmeandist = 0, gemmeandist = 0, scitmeandist = 0;
2548 int mvdcount = 0, sttcount = 0, gemcount = 0, scitcount = 0;
2549 for(
int ihit = 0; ihit < skewhitlist2.
GetNofHits(); ihit++) {
2554 mvdmeandist += zdistance;
2557 else if(hit->
IsStt()) {
2558 sttmeandist += zdistance;
2561 else if(hit->
IsGem()) {
2562 gemmeandist += zdistance;
2566 scitmeandist += zdistance;
2570 if(mvdcount > 0) mvdmeandist /= mvdcount;
2571 if(sttcount > 0) sttmeandist /= sttcount;
2572 if(gemcount > 0) gemmeandist /= gemcount;
2573 if(scitcount > 0) scitmeandist /= scitcount;
2576 for(
int ihit = 0; ihit < skewhitlist2.
GetNofHits(); ihit++) {
2582 if(
fabs(mvdmeandist - zdistance) < 2) skewhitlist3.
AddHit(hit);
2584 else if(hit->
IsStt()) {
2586 if(
fabs(sttmeandist - zdistance) < 10) skewhitlist3.
AddHit(hit);
2588 else if(hit->
IsGem()) {
2590 if(
fabs(gemmeandist - zdistance) < 20) skewhitlist3.
AddHit(hit);
2594 if(
fabs(scitmeandist - zdistance) < 20) skewhitlist3.
AddHit(hit);
2605 cout <<
"final z " << skewhitlist3.
GetNofHits() << endl;
2606 for(
int ihit = 0; ihit < skewhitlist3.
GetNofHits(); ihit++) {
2612 mrkzphi->SetMarkerColor(kMagenta);
2613 mrkzphi->Draw(
"SAME");
2626 for(
int ihit = 0; ihit < skewhitlist3.
GetNofHits(); ihit++) {
2627 hit = skewhitlist3.
GetHit(ihit);
2635 TMarker *mrkzphi =
new TMarker(phi, position.Z(), 20);
2636 mrkzphi->SetMarkerColor(kBlue);
2637 if(hit->
IsStt() == kTRUE) mrkzphi->SetMarkerColor(kYellow);
2638 mrkzphi->Draw(
"SAME");
2644 double fitm4, fitq4;
2649 TLine *line =
new TLine(0, fitq4, 360, 360 * fitm4 + fitq4);
2650 line->SetLineColor(4);
2659 for(
int ihit = 0; ihit < skewhitlist2.
GetNofHits(); ihit++) {
2660 hit = skewhitlist2.
GetHit(ihit);
2689 finaltrack.
SetZ0(z0);
2729 std::vector< int > merged;
2732 std::map < int , std::vector < int > > mergingtracks;
2735 if(merged[itrk] == 1)
continue;
2738 std::vector< int > mtracks;
2740 if(merged[jtrk] == 1)
continue;
2743 if(clusterj.
SharedAt(&clusteri, 0.70) == kTRUE) {
2744 mtracks.push_back(jtrk);
2750 if(mtracks.size() == 0) mtracks.push_back(-1);
2751 mergingtracks.insert(std::pair<
int , std::vector < int > >(itrk, mtracks));
2757 std::map< int, std::vector < int > >::iterator it;
2761 it = mergingtracks.find(itrk);
2762 if(it == mergingtracks.end())
continue;
2764 std::vector< int > mtracks = it->second;
2766 if(mtracks[0] != -1) {
2769 for(
size_t ktrk = 0; ktrk < mtracks.size(); ktrk++) {
2770 int jtrk = mtracks[ktrk];
2781 double delta, trasl[2];
2793 double rc_of_min, rc_of_max;
2797 double u = chit->
GetU();
2798 double v = chit->
GetV();
2813 double rimin, rimax;
2814 r1 < r2 ? (rimin =
r1, rimax =
r2) : (rimin = r2, rimax = r1);
2816 rimin <
fRmin ? (rc_of_min = rc,
fRmin = rimin) : fRmin;
2817 rimax >
fRmax ? (rc_of_max = rc,
fRmax = rimax) : fRmax;
2826 double delt = 0.5 *
fabs(dv - du);
2827 du < dv ? (fUmin -= delt,
fUmax += delt) : (fVmin -= delt,
fVmax += delt);
2832 for(
int ihit = 0; ihit < nofconfhits; ihit++) {
2847 double xc3, yc3, R3;
2850 bool fitting =
AnalyticalFit(&clusteri, xc2, yc2, R2, fitm, fitq);
2851 if(fitting == kFALSE)
continue;
2856 TLine *line =
new TLine(-10, -10 * fitm + fitq, 10, 10 * fitm + fitq);
2857 line->SetLineColor(kMagenta);
2862 double fita, fitb, fitc, epsilon;
2864 if(fitting == kFALSE)
continue;
2874 TArc *arcm =
new TArc(xc3, yc3, R3);
2875 arcm->SetFillStyle(0);
2876 arcm->SetLineColor(kMagenta);
2896 for(
int jtrk = 0; jtrk < mergedtracklist->
GetNofTracks(); jtrk++) {
2915 for(
int itrk = 0; itrk < mergedtracklist->
GetNofTracks(); itrk++) {
2923 double gap[5] = {0, 0, 0, 0, 0};
2924 int gapcounter[5] = {0, 0, 0, 0, 0};
2928 hit = clusteri.
GetHit(0);
2943 for(
int ihit = 0; ihit < clusteri.
GetNofHits(); ihit++) {
2944 hit = clusteri.
GetHit(ihit);
2949 if(detID == FairRootManager::Instance()->GetBranchId(
fMvdPixelBranch)) {
2950 gap[0] +=
fabs(
TMath::Sqrt((pos.X() - xc3) * (pos.X() - xc3) + (pos.Y() - yc3) * (pos.Y() - yc3)) - R3);
2954 else if(detID == FairRootManager::Instance()->GetBranchId(
fMvdStripBranch)) {
2955 gap[1] +=
fabs(
TMath::Sqrt((pos.X() - xc3) * (pos.X() - xc3) + (pos.Y() - yc3) * (pos.Y() - yc3)) - R3);
2959 else if(detID == FairRootManager::Instance()->GetBranchId(
fSttBranch)) {
2964 else if(detID == FairRootManager::Instance()->GetBranchId(
fSciTBranch)) {
2965 gap[3] +=
fabs(
TMath::Sqrt((pos.X() - xc3) * (pos.X() - xc3) + (pos.Y() - yc3) * (pos.Y() - yc3)) - R3);
2969 else if(detID == FairRootManager::Instance()->GetBranchId(
fGemBranch)) {
2970 gap[4] +=
fabs(
TMath::Sqrt((pos.X() - xc3) * (pos.X() - xc3) + (pos.Y() - yc3) * (pos.Y() - yc3)) - R3);
2976 for(
int igap = 0; igap < 5; igap++) {
2977 if(gapcounter[igap] > 0) {
2978 gap[igap]/= gapcounter[igap];
2979 gap[igap] =
fabs(gap[igap]);
2982 if(gap[0] > 2) vote =
false;
2983 if(gap[1] > 2) vote =
false;
2984 if(gap[2] > 10) vote =
false;
2985 if(gap[3] > 2) vote =
false;
2986 if(gap[4] > 2) vote =
false;
2989 if(vote ==
true) cleanedtracklist->
AddTrack(track);
2991 delete mergedtracklist;
2996 for(
int itrk = 0; itrk < cleanedtracklist->
GetNofTracks(); itrk++) {
3003 clusteri.
Draw(kRed);
3008 for(
int ihit = 0; ihit < clusteri.
GetNofHits(); ihit++) {
3009 hit = clusteri.
GetHit(ihit);
3018 for(
int jhit = 0; jhit < indiv2.GetEntriesFast(); jhit++) {
3045 clusteri.
Draw(kMagenta);
3053 delete cleanedtracklist;
3055 for(
int itrk = 0; itrk < indivtracklist->
GetNofTracks(); itrk++) {
3059 delete indivtracklist;
3080 for(
int ihit = 0; ihit < tempcluster.
GetNofHits(); ihit++) {
3081 hit = tempcluster.
GetHit(ihit);
3082 if(hit->
IsSciTil() == kTRUE)
continue;
3083 tempcluster2.
AddHit(hit);
3091 Int_t size = clref1.GetEntriesFast();
3097 size = clref2.GetEntriesFast();
3101 size = clref3.GetEntriesFast();
3117 cout <<
"TRACK " << itrk << endl;
3118 cout <<
"MOM FIRST: TOT, PT, PL " << outputtrack->
GetParamFirst().GetMomentum().Mag() <<
" " << outputtrack->
GetParamFirst().GetMomentum().Perp() <<
" " << outputtrack->
GetParamFirst().GetMomentum().Z() <<
" nofhits " << outputtrackcand->
GetNHits() << endl;
3120 cout <<
"Z0, TANL " << track->
GetZ0() <<
" " << track->
GetTanL() << endl;
3121 cout <<
"CHARGE " << track->
GetCharge() << endl;
3127 int noflongtracks =
fTrackArray->GetEntriesFast();
3144 cout <<
"ERROR track " << itrk <<
" has no candidate association" << endl;
3148 for (
size_t ihit = 0; ihit < cand->
GetNHits(); ihit++) {
3153 if(detId == FairRootManager::Instance()->GetBranchId(
fMvdPixelBranch)) {
3157 else if(detId == FairRootManager::Instance()->GetBranchId(
fMvdStripBranch)) {
3161 else if(detId == FairRootManager::Instance()->GetBranchId(
fSttBranch)) {
3165 else if(detId == FairRootManager::Instance()->GetBranchId(
fSciTBranch)) {
3169 else if(detId == FairRootManager::Instance()->GetBranchId(
fGemBranch)) {
3180 for(
int itrk = 0; itrk <
fTrackArray->GetEntriesFast(); itrk++) {
3198 int nofgemhits_on_layer[
NOFLAYERS] = {0, 0, 0, 0, 0, 0};
3201 if(hit->
IsUsed() == kTRUE)
continue;
3205 nofgemhits_on_layer[layerid]++;
3209 for(
int ihit = 0; ihit < gemcluster.
GetNofHits(); ihit++) {
3218 for(
int ihit = 0; ihit < gemcluster.
GetNofHits(); ihit++) {
3220 if(hiti->
IsUsed() == kTRUE)
continue;
3240 while(goOn ==
true) {
3242 double tmpdistance = 1000;
3243 for(
int jhit = 0; jhit < gemcluster.
GetNofHits(); jhit++) {
3245 if(hitj->
IsUsed() == kTRUE)
continue;
3248 if(layerjd != tmphitj->
GetSensorID() + 1)
continue;
3251 if(distance < tmpdistance) {
3253 tmpdistance = distance;
3258 if(tmpdistance < 20) {
3265 tmphitj->
Draw(kBlue);
3320 Double_t delta = 0, trasl[2] = {0., 0.};
3323 int tmpsensid = 1000;
3324 for(
int jhit = 0; jhit < cluster.
GetNofHits(); jhit++) {
3342 double rc_of_min, rc_of_max;
3346 double u = chit->
GetU();
3347 double v = chit->
GetV();
3351 double sigma = 1000.;
3353 if(TMath::IsNaN(chit->
GetPosition().X()))
continue;
3368 double rimin, rimax;
3369 r1 < r2 ? (rimin =
r1, rimax =
r2) : (rimin = r2, rimax = r1);
3371 rimin <
fRmin ? (rc_of_min = rc,
fRmin = rimin) : fRmin;
3372 rimax >
fRmax ? (rc_of_max = rc,
fRmax = rimax) : fRmax;
3381 double delt =
fabs(dv - du)/2.;
3382 du < dv ? (fUmin -= delt,
fUmax += delt) : (fVmin -= delt,
fVmax += delt);
3387 for(
int khit = 0; khit < nofconfhits; khit++) {
3412 TArc *arcm =
new TArc(xc, yc, R);
3413 arcm->SetFillStyle(0);
3414 arcm->SetLineColor(3);
3425 std::vector< std::pair< int, int > > dontuse;
3429 std::map< int, int > sectorids;
3430 for(
int jhit = 0; jhit < cluster.
GetNofHits(); jhit++)
3432 hit = cluster.
GetHit(jhit);
3436 if(sectorids.find(secid) == sectorids.end()) sectorids[secid] = 1;
3437 else sectorids[secid]++;
3440 int tmpsecentries = 0, tmpsec = -1;
3441 for(
size_t isec = 0;
isec < sectorids.size();
isec++) {
3442 if(tmpsecentries < sectorids[
isec]) {
3443 tmpsecentries = sectorids[
isec];
3447 int sectorID = tmpsec;
3451 bool border =
false;
3452 int othersecID = -1;
3453 for(
size_t isec = 0;
isec < sectorids.size();
isec++) {
3454 if(sectorids[
isec] > 0 && (
int)
isec != sectorID) {
3465 double xydistance =
TMath::Sqrt(xc * xc + yc * yc);
3468 bool searchnearby =
false;
3475 searchnearby =
true;
3482 searchnearby =
true;
3486 searchnearby =
true;
3494 bool crossL =
false;
3495 double y1left = 0, y2left = 0;
3502 else if(delta1 > 0) {
3508 bool crossR =
false;
3509 double y1right = 0, y2right = 0;
3515 else if(delta2 > 0) {
3521 bool crossingpipe =
false;
3522 if(crossL && crossR) {
3523 double ycross = 0.5 * (y1left + y1right);
3525 ycross = 0.5 * (y2left + y2right);
3534 crossingpipe =
true;
3546 if(hit->
IsUsed() == kTRUE)
continue;
3549 if(crossingpipe ==
true) {
3550 if(((hit->
GetSector() == 0 && sectorID == 5) && !(hit->
GetSector() == 2 && sectorID == 3) && !(hit->
GetSector() == 5 && sectorID == 0) && !(hit->
GetSector() == 3 && sectorID == 2)) && searchnearby ==
false)
continue;
3552 else if(searchnearby ==
true) {
3556 if(border ==
false && hit->
GetSector() != sectorID)
continue;
3557 else if(border ==
true && (hit->
GetSector() != sectorID && hit->
GetSector() != othersecID))
continue;
3561 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
3562 double recoiso =
fabs(distance_hit_center - R);
3564 if(recoiso < 3.) cluster2.
AddHit(hit);
3571 for(
int jhit = 0; jhit < cluster2.
GetNofHits(); jhit++)
3573 hit = cluster2.
GetHit(jhit);
3577 if(sectorids.find(secid) == sectorids.end()) sectorids[secid] = 1;
3578 else sectorids[secid]++;
3581 tmpsecentries = 0, tmpsec = -1;
3582 for(
size_t isec = 0;
isec < sectorids.size();
isec++) {
3583 if(tmpsecentries < sectorids[
isec]) {
3584 tmpsecentries = sectorids[
isec];
3594 for(
size_t isec = 0;
isec < sectorids.size();
isec++) {
3595 if(sectorids[
isec] > 0 && (
int)
isec != sectorID) {
3607 if(hit->
IsUsed() == kTRUE)
continue;
3608 if(searchnearby ==
true) {
3612 if(border ==
false && hit->
GetSector() != sectorID)
continue;
3613 else if(border ==
true && (hit->
GetSector() != sectorID && hit->
GetSector() != othersecID))
continue;
3616 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
3618 double recoiso =
fabs(distance_hit_center - R);
3623 int samesensor =
false;
3624 for(
int khit = 0; khit < cluster2.
GetNofHits(); khit++) {
3629 double distance_hitk_center = (hitk->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
3630 double recoisok =
fabs(distance_hitk_center - R);
3631 if(recoiso < recoisok) {
3634 dontuse.push_back(dontpair);
3639 dontuse.push_back(dontpair);
3643 if(samesensor ==
false) cluster2.
AddHit(hit);
3658 if(hit->
IsUsed() == kTRUE)
continue;
3659 if(searchnearby ==
true) {
3663 if(border ==
false && hit->
GetSector() != sectorID)
continue;
3664 else if(border ==
true && (hit->
GetSector() != sectorID && hit->
GetSector() != othersecID))
continue;
3667 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
3668 double recoiso =
fabs(distance_hit_center - R);
3673 int samesensor =
false;
3674 for(
int khit = 0; khit < cluster2.
GetNofHits(); khit++) {
3679 double distance_hitk_center = (hitk->
GetPosition().XYvector() - TVector2(xc, yc)).Mod();
3680 double recoisok =
fabs(distance_hitk_center - R);
3681 if(recoiso < recoisok) {
3684 dontuse.push_back(dontpair);
3689 dontuse.push_back(dontpair);
3693 if(samesensor ==
false) cluster2.
AddHit(hit);
3708 cluster2.
Draw(kRed);
3712 cout <<
"want to go on?" << endl;
3720 for(
int jhit = 0; jhit < cluster2.
GetNofHits(); jhit++) {
3721 hit = cluster2.
GetHit(jhit);
3728 for(
int jhit = 0; jhit < cluster2.
GetNofHits(); jhit++) {
3729 hit = cluster2.
GetHit(jhit);
3756 for(
int jhit = 0; jhit < cluster2.
GetNofHits(); jhit++) {
3757 hit = cluster2.
GetHit(jhit);
3766 for(
int jhit = 0; jhit < cluster2.
GetNofHits(); jhit++) {
3767 hit = cluster2.
GetHit(jhit);
3768 if(hit->
IsStt())
continue;
3774 TMarker *mrkz = NULL;
3775 if(hit->
IsMvdPixel()) mrkz =
new TMarker(phi, position.Z(), 21);
3776 else if(hit->
IsMvdStrip()) mrkz =
new TMarker(phi, position.Z(), 25);
3777 else if(hit->
IsStt()) mrkz =
new TMarker(phi, position.Z(), 6);
3778 else if(hit->
IsGem()) mrkz =
new TMarker(phi, position.Z(), 24);
3779 else if(hit->
IsSciTil()) mrkz =
new TMarker(phi, position.Z(), 26);
3781 mrkz->SetMarkerColor(kBlue);
3787 for(
int khit = jhit + 1; khit < cluster2.
GetNofHits(); khit++) {
3789 if(hitk->
IsStt())
continue;
3791 double phik = hitk->
GetPhi();
3793 double denom = (phi - phik) * (phi - phik) + (positionk.Z() - position.Z()) * (positionk.Z() - position.Z());
3794 if(denom == 0)
continue;
3795 double cost = (positionk.Z() - position.Z()) /
TMath::Sqrt(denom);
3796 double theta = TMath::ACos(cost);
3797 double r1 = phi * cost + position.Z() *
TMath::Sin(theta);
3798 double r2 = phik * cost + positionk.Z() *
TMath::Sin(theta);
3800 if(
fabs(r1 - r2) > 1.e-9) {
3801 theta = -TMath::ACos(cost);
3802 r = phi * cost + position.Z() *
TMath::Sin(theta);
3804 fLineHisto->Fill(theta * TMath::RadToDeg(), r);
3809 int binx, biny, binz;
3810 fLineHisto->GetBinXYZ(bin, binx, biny, binz);
3811 double tpeak =
fLineHisto->GetXaxis()->GetBinCenter(binx);
3812 double rpeak =
fLineHisto->GetYaxis()->GetBinCenter(biny);
3817 for(
int jhit = 0; jhit < cluster2.
GetNofHits(); jhit++) {
3818 hit = cluster2.
GetHit(jhit);
3819 if(!hit->
IsGem())
continue;
3826 TMarker *mrkz = NULL;
3827 if(hit->
IsMvdPixel()) mrkz =
new TMarker(phi, position.Z(), 21);
3828 else if(hit->
IsMvdStrip()) mrkz =
new TMarker(phi, position.Z(), 25);
3829 else if(hit->
IsStt()) mrkz =
new TMarker(phi, position.Z(), 6);
3830 else if(hit->
IsGem()) mrkz =
new TMarker(phi, position.Z(), 24);
3831 else if(hit->
IsSciTil()) mrkz =
new TMarker(phi, position.Z(), 26);
3833 mrkz->SetMarkerColor(kBlue);
3840 double sigma = 1e-5 * position.Z();
3845 double fitm4, fitq4;
3852 double fitq44 = rpeak/
TMath::Sin(tpeak * TMath::DegToRad());
3858 TLine *line =
new TLine(-360, -360 * fitm4 + fitq4, 360, 360 * fitm4 + fitq4);
3859 line->SetLineColor(3);
3862 TLine *line4 =
new TLine(-360, -360 * fitm44 + fitq44, 360, 360 * fitm44 + fitq44);
3863 line4->SetLineColor(4);
3864 line4->Draw(
"SAME");
3875 for(
int jhit = 0; jhit < cluster2.
GetNofHits(); jhit++) {
3876 hit = cluster2.
GetHit(jhit);
3877 if(hit->
IsStt())
continue;
3880 double distance =
fabs(position.Z() - fitm4 * phi - fitq4);
3883 if(distance > 100) {
3885 double comp_phi = (position.Z() - fitq4)/fitm4;
3886 if(
fabs(comp_phi - phi) < 361 &&
fabs(comp_phi - phi) > 359){
3887 if(comp_phi > phi) phi += 360.;
3890 distance =
fabs(position.Z() - fitm4 * phi - fitq4);
3899 for(
int jhit = 0; jhit < cluster2.
GetNofHits(); jhit++) {
3900 hit = cluster2.
GetHit(jhit);
3910 double distance =
fabs(position.Z() - fitm4 * phi - fitq4);
3913 if(distance > 3 && hit->
IsMvd() == kTRUE) {
3919 std::pair< int, int > dontpair(hitId, hit->
GetDetectorID());
3920 dontuse.push_back(dontpair);
3926 else if(distance > 20 && hit->
IsGem() == kTRUE) {
3931 std::pair< int, int > dontpair(hitId, hit->
GetDetectorID());
3932 dontuse.push_back(dontpair);
3942 TMarker *mrkz =
new TMarker(phi, position.Z(), 20);
3943 mrkz->SetMarkerColor(kOrange);
3953 bool stat[3] = {
false,
false,
false};
3954 for(
int jhit = 0; jhit < cluster3.
GetNofHits(); jhit++) {
3955 hit = cluster3.
GetHit(jhit);
3956 if(hit->
IsGem() == kFALSE)
continue;
3973 if(stat[0] ==
false && stat[1] ==
false && stat[2] ==
false)
continue;
3980 for(
int jhit = 0; jhit < cluster3.
GetNofHits(); jhit++) {
3981 hit = cluster3.
GetHit(jhit);
3982 if(hit->
IsStt())
continue;
3989 TMarker *mrkz = NULL;
3990 if(hit->
IsMvdPixel()) mrkz =
new TMarker(phi, position.Z(), 21);
3991 else if(hit->
IsMvdStrip()) mrkz =
new TMarker(phi, position.Z(), 25);
3992 else if(hit->
IsStt()) mrkz =
new TMarker(phi, position.Z(), 6);
3993 else if(hit->
IsGem()) mrkz =
new TMarker(phi, position.Z(), 24);
3995 mrkz->SetMarkerColor(kBlue);
4005 double fitm5, fitq5;
4010 TLine *line =
new TLine(-360, -360 * fitm5 + fitq5, 360, 360 * fitm5 + fitq5);
4011 line->SetLineColor(4);
4026 refhit = cluster3.
GetHit(0);
4042 fUmin = 1000, fVmin = 1000,
fRmin = 1000;
4047 double u = chit->
GetU();
4048 double v = chit->
GetV();
4050 if(TMath::IsNaN(u))
continue;
4053 u - rc < fUmin ? fUmin = u - rc :
fUmin;
4054 v - rc < fVmin ? fVmin = v - rc :
fVmin;
4055 u + rc > fUmax ? fUmax = u + rc :
fUmax;
4064 double rimin, rimax;
4065 r1 < r2 ? (rimin =
r1, rimax =
r2) : (rimin = r2, rimax = r1);
4067 rimin <
fRmin ? (rc_of_min = rc,
fRmin = rimin) : fRmin;
4068 rimax >
fRmax ? (rc_of_max = rc,
fRmax = rimax) : fRmax;
4077 delt =
fabs(dv - du)/2.;
4078 du < dv ? (fUmin -= delt, fUmax += delt) : (fVmin -= delt,
fVmax += delt);
4083 for(
int jhit = 0; jhit < nofconfhits; jhit++) {
4090 double fitm3, fitq3;
4091 bool fitting =
AnalyticalFit(&cluster3, xc, yc, R, fitm3, fitq3);
4092 if(fitting == kFALSE)
continue;
4093 double xc3, yc3, R3;
4099 TArc *arcm =
new TArc(xc3, yc3, R3);
4100 arcm->SetFillStyle(0);
4101 arcm->SetLineColor(4);
4119 std::vector < std::pair< int, int > >::iterator itr = dontuse.begin();
4120 itr = std::find(dontuse.begin(), dontuse.end(), thispair);
4121 if(itr != dontuse.end())
continue;
4123 if(border ==
false && hit->
GetSector() != sectorID)
continue;
4124 else if(border ==
true && (hit->
GetSector() != sectorID && hit->
GetSector() != othersecID))
continue;
4125 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc3, yc3)).Mod();
4126 double recoiso =
fabs(distance_hit_center - R3);
4131 int samesensor =
false;
4134 if(!hitk->
IsGem())
continue;
4137 double distance_hitk_center = (hitk->
GetPosition().XYvector() - TVector2(xc3, yc3)).Mod();
4138 double recoisok =
fabs(distance_hitk_center - R3);
4139 if(recoiso < recoisok) {
4142 dontuse.push_back(dontpair);
4147 dontuse.push_back(dontpair);
4170 int laststat = -1, almostlaststat = -1;
4174 if(laststat != -1 && almostlaststat == -1 && hit->
GetSensorID() != laststat) almostlaststat = hit->
GetSensorID();
4177 int hitsonlaststat = 0;
4187 if(laststat == -1 || almostlaststat == -1)
continue;
4190 if(hitsonlaststat > 1) {
4191 double dist45 = 1000;
4198 if(hitk->
GetSensorID() != almostlaststat)
continue;
4199 else if(hit->
GetSensorID() < almostlaststat)
break;
4201 if(tmpdist < dist45) fromhere = hit->
GetPosition();
4217 if(border ==
false && tube->
GetSectorID() != sectorID)
continue;
4218 else if(border ==
true && (tube->
GetSectorID() != sectorID && tube->
GetSectorID() != othersecID))
continue;
4219 double distance_hit_center = (tube->
GetPosition().XYvector() - TVector2(xc3, yc3)).Mod();
4220 double recoiso =
fabs(distance_hit_center - R3);
4234 std::vector < std::pair< int, int > >::iterator itr = dontuse.begin();
4235 itr = std::find(dontuse.begin(), dontuse.end(), thispair);
4236 if(itr != dontuse.end())
continue;
4238 if(border ==
false && hit->
GetSector() != sectorID)
continue;
4239 else if(border ==
true && (hit->
GetSector() != sectorID && hit->
GetSector() != othersecID))
continue;
4240 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc3, yc3)).Mod();
4241 double recoiso =
fabs(distance_hit_center - R3);
4246 int samesensor =
false;
4252 double distance_hitk_center = (hitk->
GetPosition().XYvector() - TVector2(xc3, yc3)).Mod();
4253 double recoisok =
fabs(distance_hitk_center - R3);
4254 if(recoiso < recoisok) {
4257 dontuse.push_back(dontpair);
4262 dontuse.push_back(dontpair);
4280 std::vector < std::pair< int, int > >::iterator itr = dontuse.begin();
4281 itr = std::find(dontuse.begin(), dontuse.end(), thispair);
4282 if(itr != dontuse.end())
continue;
4284 if(border ==
false && hit->
GetSector() != sectorID)
continue;
4285 else if(border ==
true && (hit->
GetSector() != sectorID && hit->
GetSector() != othersecID))
continue;
4286 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc3, yc3)).Mod();
4287 double recoiso =
fabs(distance_hit_center - R3);
4298 double tmpdistance = 1000;
4301 double distance_hit_center = (hit->
GetPosition().XYvector() - TVector2(xc3, yc3)).Mod();
4302 double recoiso =
fabs(distance_hit_center - R3);
4304 if(recoiso < tmpdistance) {
4306 tmpdistance = recoiso;
4328 if(phi == 360) phi = 0;
4347 mrk->SetMarkerColor(kMagenta);
4360 if(!hit->
IsGem())
continue;
4371 if(!hit->
IsGem())
continue;
4373 if(phi < 10) phi += 360;
4386 mrk->SetMarkerColor(kGreen);
4407 mrk->SetMarkerColor(kOrange);
4433 fUmin = 1000, fVmin = 1000,
fRmin = 1000;
4438 double u = chit->
GetU();
4439 double v = chit->
GetV();
4441 if(TMath::IsNaN(u))
continue;
4444 u - rc < fUmin ? fUmin = u - rc :
fUmin;
4445 v - rc < fVmin ? fVmin = v - rc :
fVmin;
4446 u + rc > fUmax ? fUmax = u + rc :
fUmax;
4455 double rimin, rimax;
4456 r1 < r2 ? (rimin =
r1, rimax =
r2) : (rimin = r2, rimax = r1);
4458 rimin <
fRmin ? (rc_of_min = rc,
fRmin = rimin) : fRmin;
4459 rimax >
fRmax ? (rc_of_max = rc,
fRmax = rimax) : fRmax;
4468 delt =
fabs(dv - du)/2.;
4469 du < dv ? (fUmin -= delt, fUmax += delt) : (fVmin -= delt,
fVmax += delt);
4474 for(
int jhit = 0; jhit < nofconfhits; jhit++) {
4482 double fitm6, fitq6;
4484 if(fitting == kFALSE)
continue;
4486 double xc6, yc6, R6;
4490 if(R6 > 2500)
continue;
4498 TLine *line =
new TLine(-10, -10 * fitm6 + fitq6, 10, 10 * fitm6 + fitq6);
4499 line->SetLineColor(kMagenta);
4503 TArc *arcm =
new TArc(xc6, yc6, R6);
4504 arcm->SetFillStyle(0);
4505 arcm->SetLineColor(kMagenta);
4537 if(hit->
IsSttSkew() == kFALSE)
continue;
4539 if(border ==
false && hit->
GetSector() != sectorID)
continue;
4540 else if(border ==
true && (hit->
GetSector() != sectorID && hit->
GetSector() != othersecID))
continue;
4548 TVector3 first = tube->
GetPosition() + wireDirection * halflength;
4549 TVector3 second = tube->
GetPosition() - wireDirection * halflength;
4555 TVector2 intersection1, intersection2;
4558 if(nofintersections == 0)
continue;
4559 if(nofintersections >= 2) {
4560 cout <<
"ERROR: MORE THAN 1 INTERSECTION!!" << endl;
4571 double beta = wireDirection.Phi();
4576 TVector2 rottangent(rtx, rty);
4577 rottangent = rottangent.Unit();
4588 Double_t rotm = rottangent.Y()/rottangent.X();
4598 x0a = (-rotp +
TMath::Sqrt(b * b + a * a * rotm * rotm)) / rotm;
4599 x0b = (-rotp -
TMath::Sqrt(b * b + a * a * rotm * rotm)) / rotm;
4602 double intxa = (x0a * b * b - rotm * rotp * a *
a) / (b * b + rotm * rotm * a * a);
4603 double intya = rotm * intxa + rotp;
4604 double intxb = (x0b * b * b - rotm * rotp * a *
a) / (b * b + rotm * rotm * a * a);
4605 double intyb = rotm * intxb + rotp;
4632 double y0a = y0anew;
4634 double y0b = y0bnew;
4640 TEllipse *ell1 =
new TEllipse(x0a, y0a, a, b, 0, 360, -beta);
4641 ell1->SetFillStyle(0);
4642 ell1->SetLineColor(4);
4644 TEllipse *ell2 =
new TEllipse(x0b, y0b, a, b, 0, 360, -beta);
4645 ell2->SetFillStyle(0);
4646 ell2->SetLineColor(6);
4649 TMarker *mrkinta =
new TMarker(intxa, intya, 20);
4650 mrkinta->SetMarkerColor(4);
4651 mrkinta->Draw(
"SAME");
4652 TMarker *mrkintb =
new TMarker(intxb, intyb, 20);
4653 mrkintb->SetMarkerColor(6);
4654 mrkintb->Draw(
"SAME");
4664 Double_t t = ((x0a + y0a) - (first.X() + first.Y())) / ((second.X() - first.X()) + (second.Y() - first.Y()));
4665 Double_t z0a = first.Z() + (second.Z() - first.Z()) * t;
4668 t = ((x0b + y0b) - (first.X() + first.Y())) / ((second.X() - first.X()) + (second.Y() - first.Y()));
4669 Double_t z0b = first.Z() + (second.Z() - first.Z()) * t;
4671 TVector3 center1(x0a, y0a, z0a);
4672 TVector3 center2(x0b, y0b, z0b);
4675 double dx = intxa - x0a;
4676 double dy = intya - y0a;
4677 TVector3 dxdy(dx, dy, 0.0);
4679 TVector3 tfirst = first + dxdy;
4680 TVector3 tsecond = second + dxdy;
4682 t = ((intxa + intya) - (tfirst.X() + tfirst.Y())) / ((tsecond.X() - tfirst.X()) + (tsecond.Y() - tfirst.Y()));
4683 double intza = tfirst.Z() + (tsecond.Z() - tfirst.Z()) * t;
4687 TLine *linezx1 =
new TLine(tfirst.X(), tfirst.Z(), tsecond.X(), tsecond.Z());
4688 linezx1->SetLineStyle(1);
4689 linezx1->Draw(
"SAME");
4690 TMarker *mrkza1 =
new TMarker(intxa, intza, 20);
4691 mrkza1->SetMarkerColor(kBlue - 9);
4692 mrkza1->Draw(
"SAME");
4696 tfirst = first - dxdy;
4697 tsecond = second - dxdy;
4699 t = ((intxb + intyb) - (tfirst.X() + tfirst.Y())) / ((tsecond.X() - tfirst.X()) + (tsecond.Y() - tfirst.Y()));
4700 double intzb = tfirst.Z() + (tsecond.Z() - tfirst.Z()) * t;
4702 TVector3 fin_intersection1(intxa, intya, intza);
4703 TVector3 fin_intersection2(intxb, intyb, intzb);
4707 double phi1 = finaltrack.
ComputePhiFrom(fin_intersection1, fromhere);
4708 double phi2 = finaltrack.
ComputePhiFrom(fin_intersection2, fromhere);
4709 hit->
SetPhi((phi1 + phi2) * 0.5);
4711 skewhitlist.
AddHit(
PndTrkSkewHit(*hit, trackID, center1, fin_intersection1, phi1, center2, fin_intersection2, phi2, a, b, -1, beta));
4720 if(hit->
IsStt() == kFALSE) {
4735 mrk->SetMarkerColor(kBlue);
4745 double fitm7 = fitm5;
4747 for(
int jhit = 0; jhit < skewhitlist.
GetNofHits(); jhit++) {
4748 hit = skewhitlist.
GetHit(jhit);
4750 if(hit->
IsStt() == kFALSE) {
4759 for(
int jhit = 0; jhit < skewhitlist.
GetNofHits(); jhit++) {
4760 hit = skewhitlist.
GetHit(jhit);
4763 if(hit->
IsStt() == kTRUE) {
4768 double phi1 = skewhit->
GetPhi1();
4769 double phi2 = skewhit->
GetPhi2();
4772 TLine *linezphi =
new TLine(phi1, fin_intersection1.Z(), phi2, fin_intersection2.Z());
4773 linezphi->SetLineStyle(1);
4774 linezphi->Draw(
"SAME");
4775 TMarker *mrkzphi1 =
new TMarker(phi1, fin_intersection1.Z(), 20);
4776 mrkzphi1->SetMarkerColor(kBlue - 9);
4777 mrkzphi1->Draw(
"SAME");
4778 TMarker *mrkzphi2 =
new TMarker(phi2, fin_intersection2.Z(), 20);
4779 mrkzphi2->SetMarkerColor(kMagenta - 7);
4780 mrkzphi2->Draw(
"SAME");
4787 mrkzphi->SetMarkerColor(kOrange);
4788 mrkzphi->Draw(
"SAME");
4799 TLine *line =
new TLine(-360, -360 * fitm7 + fitq7, 360, 360 * fitm7 + fitq7);
4800 line->SetLineColor(2);
4811 for(
int jhit = 0; jhit < skewhitlist.
GetNofHits(); jhit++) {
4815 skewhitlist2.
AddHit(hit);
4819 TVector3 fin_intersection1 = ((
PndTrkSkewHit*) hit)->GetIntersection1();
4820 TVector3 fin_intersection2 = ((
PndTrkSkewHit*) hit)->GetIntersection2();
4824 double dist1 =
fabs(fitm7 * phi1 - fin_intersection1.Z() + fitq7) /
TMath::Sqrt(fitm7 * fitm7 + 1);
4825 double dist2 =
fabs(fitm7 * phi2 - fin_intersection2.Z() + fitq7) /
TMath::Sqrt(fitm7 * fitm7 + 1);
4827 double distance = 1000;
4828 dist1 < dist2 ? (distance = dist1, hit->
SetPosition(fin_intersection1), hit->
SetPhi(phi1)) : (distance = dist2, hit->
SetPosition(fin_intersection2), hit->
SetPhi(phi2));
4830 if(distance < 10) skewhitlist2.
AddHit(hit);
4835 TVector3 position, positionk;
4837 for(
int jhit = 0; jhit < skewhitlist2.
GetNofHits(); jhit++) {
4838 hit = skewhitlist2.
GetHit(jhit);
4842 if(hit->
IsStt()) solutions = 2;
4845 for(
int khit = jhit + 1; khit < skewhitlist2.
GetNofHits(); khit++) {
4850 if(hitk->
IsStt()) ksolutions = 2;
4851 else ksolutions = 1;
4853 for(
int jsol = 0; jsol < solutions; jsol++)
4855 if(solutions == 1) {
4859 else if(jsol == 0) {
4868 for(
int ksol = 0; ksol < ksolutions; ksol++)
4870 if(ksolutions == 1) {
4874 else if(ksol == 0) {
4883 double cost = (positionk.Z() - position.Z()) /
TMath::Sqrt((phi - phik) * (phi - phik) + (positionk.Z() - position.Z()) * (positionk.Z() - position.Z()));
4884 double theta = TMath::ACos(cost);
4885 double r1 = phi * cost + position.Z() *
TMath::Sin(theta);
4886 double r2 = phik * cost + positionk.Z() *
TMath::Sin(theta);
4888 if(
fabs(r1 - r2) > 1.e-9) {
4889 theta = -TMath::ACos(cost);
4890 r = phi * cost + position.Z() *
TMath::Sin(theta);
4892 fLineHisto->Fill(theta * TMath::RadToDeg(), r);
4899 fLineHisto->GetBinXYZ(bin, binx, biny, binz);
4900 tpeak =
fLineHisto->GetXaxis()->GetBinCenter(binx);
4901 rpeak =
fLineHisto->GetYaxis()->GetBinCenter(biny);
4904 double fitq8 = rpeak/
TMath::Sin(tpeak * TMath::DegToRad());
4910 TLine *line =
new TLine(-360, -360 * fitm8 + fitq8, 360, 360 * fitm8 + fitq8);
4911 line->SetLineColor(3);
4920 for(
int jhit = 0; jhit < skewhitlist2.
GetNofHits(); jhit++) {
4921 hit = skewhitlist2.
GetHit(jhit);
4944 finalfwdtrack.
SetZ0(z0);
4967 Int_t size = clref1.GetEntriesFast();
4971 size = clref2.GetEntriesFast();
4975 size = clref3.GetEntriesFast();
4987 track->
Draw(kGreen);
4990 cout <<
"TRACK " << itrk << endl;
4991 cout <<
"MOM FIRST: TOT, PT, PL " << outputtrack->
GetParamFirst().GetMomentum().Mag() <<
" " << outputtrack->
GetParamFirst().GetMomentum().Perp() <<
" " << outputtrack->
GetParamFirst().GetMomentum().Z() <<
" nofhits " << outputtrackcand->
GetNHits() << endl;
4996 cout <<
"CHARGE " << track->
GetCharge() << endl;
5039 cout <<
"FINISH" << endl;
5073 TH2F *
h2 =
new TH2F(
"h2",
"", 100, -70, 70, 100, -70, 70);
5082 double circle[ncircles][3];
5084 for(
int i = 0;
i < ncircles;
i++)
5099 TArc *arc =
new TArc(xc, yc, rc);
5100 arc->SetFillStyle(0);
5108 for(
int a = 0;
a < ncircles;
a++) {
5109 for(
int b =
a + 1;
b < ncircles;
b++) {
5110 for(
int c =
b + 1;
c < ncircles;
c++) {
5111 for(
int csign1 = -1; csign1 <= 1; csign1 += 2) {
5112 for(
int csign2 = -1; csign2 <= 1; csign2 += 2) {
5113 double avar = 2 * (circle[
a][0] - circle[
b][0]);
5114 double bvar = 2 * (circle[
a][1] - circle[
b][1]);
5115 double cvar = 2 * (csign1 * circle[
a][2] + csign2 * circle[
b][2]);
5116 double dvar = (circle[
a][0] * circle[
a][0] + circle[
a][1] * circle[
a][1] - circle[
a][2] * circle[
a][2]) - (circle[
b][0] * circle[
b][0] + circle[
b][1] * circle[
b][1] - circle[
b][2] * circle[
b][2]);
5118 for(
int csignp = -1; csignp <= 1; csignp += 2) {
5119 double avarp = 2 * (circle[
a][0] - circle[
c][0]);
5120 double bvarp = 2 * (circle[
a][1] - circle[
c][1]);
5121 double cvarp = 2 * (csign1 * circle[
a][2] + csignp * circle[
c][2]);
5122 double dvarp = (circle[
a][0] * circle[
a][0] + circle[
a][1] * circle[
a][1] - circle[
a][2] * circle[
a][2]) - (circle[
c][0] * circle[
c][0] + circle[
c][1] * circle[
c][1] - circle[
c][2] * circle[
c][2]);
5125 double thisq = (dvar * bvarp - dvarp * bvar) / (avar * bvarp - avarp * bvar);
5126 double thism = (bvar * cvarp - bvarp * cvar) / (avar * bvarp - avarp * bvar);
5128 double thisp = (avar * dvarp - avarp * dvar) / (avar * bvarp - avarp * bvar);
5129 double thiss = (cvar * avarp - cvarp * avar) / (avar * bvarp - avarp * bvar);
5134 double thisa = thism * thism + thiss * thiss - 1;
5135 double thisb = thism * (thisq - circle[
a][0]) + thiss * (thisp - circle[
a][1]) - csign1 * circle[
a][2];
5136 double thisc = (thisq - circle[
a][0]) * (thisq - circle[a][0]) + (thisp - circle[
a][1]) * (thisp - circle[a][1]) - circle[
a][2] * circle[
a][2];
5138 if((thisb * thisb - thisa * thisc) < 0)
continue;
5139 double r = (- thisb +
TMath::Sqrt(thisb * thisb - thisa * thisc))/ thisa;
5142 r = (- thisb -
TMath::Sqrt(thisb * thisb - thisa * thisc))/ thisa;
5153 X.push_back(thism * R[counter] + thisq);
5154 Y.push_back(thiss * R[counter] + thisp);
5184 double circle[ncircles][3];
5205 if((circle[1][1] - circle[0][1])/(circle[1][0] - circle[0][0]) - (circle[2][1] - circle[0][1])/(circle[2][0] - circle[0][0]) < 0.0001){
5206 circle[0][0] += 0.01;
5207 circle[0][1] += 0.01;
5210 for(
int a = 0;
a < ncircles;
a++) {
5211 for(
int b =
a + 1;
b < ncircles;
b++) {
5212 for(
int c =
b + 1;
c < ncircles;
c++) {
5213 Y = 0.5 * ((circle[
a][0] - circle[
b][0]) * (circle[
a][0] * circle[
a][0] - circle[
c][0] * circle[
c][0] + circle[
a][1] * circle[
a][1] - circle[
c][1] * circle[
c][1]) - (circle[
a][0] - circle[
c][0]) * (circle[
a][0] * circle[
a][0] - circle[
b][0] * circle[
b][0] + circle[
a][1] * circle[
a][1] - circle[
b][1] * circle[
b][1])) / ((circle[
a][1] - circle[c][1]) * (circle[
a][0] - circle[
b][0]) - (circle[
a][1] - circle[b][1]) * (circle[
a][0] - circle[
c][0]));
5215 if(circle[
a][0] != circle[b][0]) {
5216 X = 0.5 * (circle[
a][0] * circle[
a][0] - circle[
b][0] * circle[
b][0] + circle[
a][1] * circle[
a][1] - circle[
b][1] * circle[
b][1] ) / (circle[
a][0] - circle[b][0]) - Y * (circle[
a][1] - circle[
b][1])/(circle[
a][0] - circle[b][0]);
5218 else if(circle[
a][0] != circle[c][0]) {
5219 X = 0.5 * (circle[
a][0] * circle[
a][0] - circle[
c][0] * circle[
c][0] + circle[
a][1] * circle[
a][1] - circle[
c][1] * circle[
c][1] ) / (circle[
a][0] - circle[c][0]) - Y * (circle[
a][1] - circle[
c][1])/(circle[
a][0] - circle[c][0]);
5222 R =
TMath::Sqrt((circle[
a][0] - X) * (circle[
a][0] - X) + (circle[
a][1] - Y) * (circle[
a][1] - Y));
5232 double circle[ncircles][3];
5247 for(
int a = 0;
a < ncircles;
a++) {
5248 for(
int b =
a + 1;
b < ncircles;
b++) {
5249 for(
int c =
b + 1;
c < ncircles;
c++) {
5250 for(
int csign1 = -1; csign1 <= 1; csign1 += 2) {
5251 for(
int csign2 = -1; csign2 <= 1; csign2 += 2) {
5252 double avar = 2 * (circle[
a][0] - circle[
b][0]);
5253 double bvar = 2 * (circle[
a][1] - circle[
b][1]);
5254 double cvar = 2 * (csign1 * circle[
a][2] + csign2 * circle[
b][2]);
5255 double dvar = (circle[
a][0] * circle[
a][0] + circle[
a][1] * circle[
a][1] - circle[
a][2] * circle[
a][2]) - (circle[
b][0] * circle[
b][0] + circle[
b][1] * circle[
b][1] - circle[
b][2] * circle[
b][2]);
5257 for(
int csignp = -1; csignp <= 1; csignp += 2) {
5258 double avarp = 2 * (circle[
a][0] - circle[
c][0]);
5259 double bvarp = 2 * (circle[
a][1] - circle[
c][1]);
5260 double cvarp = 2 * (csign1 * circle[
a][2] + csignp * circle[
c][2]);
5261 double dvarp = (circle[
a][0] * circle[
a][0] + circle[
a][1] * circle[
a][1] - circle[
a][2] * circle[
a][2]) - (circle[
c][0] * circle[
c][0] + circle[
c][1] * circle[
c][1] - circle[
c][2] * circle[
c][2]);
5264 double thisq = (dvar * bvarp - dvarp * bvar) / (avar * bvarp - avarp * bvar);
5265 double thism = (bvar * cvarp - bvarp * cvar) / (avar * bvarp - avarp * bvar);
5267 double thisp = (avar * dvarp - avarp * dvar) / (avar * bvarp - avarp * bvar);
5268 double thiss = (cvar * avarp - cvarp * avar) / (avar * bvarp - avarp * bvar);
5273 double thisa = thism * thism + thiss * thiss - 1;
5274 double thisb = thism * (thisq - circle[
a][0]) + thiss * (thisp - circle[
a][1]) - csign1 * circle[
a][2];
5275 double thisc = (thisq - circle[
a][0]) * (thisq - circle[a][0]) + (thisp - circle[
a][1]) * (thisp - circle[a][1]) - circle[
a][2] * circle[
a][2];
5277 if((thisb * thisb - thisa * thisc) < 0)
continue;
5278 double r = (- thisb +
TMath::Sqrt(thisb * thisb - thisa * thisc))/ thisa;
5281 r = (- thisb -
TMath::Sqrt(thisb * thisb - thisa * thisc))/ thisa;
5292 X.push_back(thism * R[counter] + thisq);
5293 Y.push_back(thiss * R[counter] + thisp);
5332 if(
fVerbose) cout <<
"Refresh stt" << endl;
5334 if(
fVerbose) cout <<
"Refresh pixel" << endl;
5336 if(
fVerbose) cout <<
"Refresh strip" << endl;
5338 if(
fVerbose) cout <<
"Refresh scit" << endl;
5340 if(
fVerbose) cout <<
"Refresh gem" << endl;
5342 if(
fVerbose) cout <<
"Refresh stop" << endl;
5370 if(
hxy == NULL)
hxy =
new TH2F(
"hxy",
"xy plane", 110, -55, 55, 110, -55, 55);
5373 hxy->SetStats(kFALSE);
5377 for(
int itube = 1; itube <
fTubeArray->GetEntriesFast(); itube++) {
5381 arc->SetFillStyle(0);
5382 arc->SetLineColor(kCyan - 10);
5391 mrk->SetMarkerColor(kCyan - 10);
5409 if(
huv == NULL)
huv =
new TH2F(
"huv",
"uv plane", 100, x1, x2, 100, y1, y2);
5412 huv->GetXaxis()->SetLimits(x1, x2);
5413 huv->GetYaxis()->SetLimits(y1, y2);
5446 for(
int j = 0; j < neighs.GetEntriesFast(); j++) {
5449 if(neighs2.GetEntriesFast() > 2) {
5452 for(
int k = 0; k < neighs2.GetEntriesFast(); k++) {
5455 cout <<
"tubes " << tubeB <<
" "<< tubeC << endl;
5459 if(counter > 2)
continue;
5485 cout <<
"new neigh hit?" << endl;
5504 for(
int i = 0;
i < neighs.GetEntriesFast();
i++) {
5519 TArc *arc =
new TArc(u, v, r);
5520 arc->SetFillStyle(0);
5524 TMarker *mrk =
new TMarker(u, v, marker);
5539 for(
int jhit = 0; jhit < cluster->
GetNofHits(); jhit++) {
5559 if(ntot == 0)
return NULL;
5565 for(
int jhit = 0; jhit < ntot; jhit++) {
5571 if(
fVerbose > 1) cout <<
"STT hit " << jhit <<
"already used " << endl;
5580 if(tmphitid == -1)
return NULL;
5583 if(
fVerbose > 1) cout <<
"STT REFERENCE HIT " << tmphitid <<
" " << refhit->
GetIsochrone() << endl;
5597 if(
fVerbose > 1) cout <<
"already used V" << endl;
5603 if(tmphitid == -1)
return NULL;
5605 if(
fVerbose > 1) cout <<
"MVD PIXEL REFERENCE HIT " << refhit->
GetHitID() << endl;
5618 if(
fVerbose > 1) cout <<
"already used V" << endl;
5624 if(tmphitid == -1)
return NULL;
5626 if(
fVerbose > 1) cout <<
"MVD STRIP REFERENCE HIT " << refhit->
GetHitID() << endl;
5635 if(refhit != NULL)
return refhit;
5646 if(refhit != NULL)
return refhit;
5657 if(ntot == 0)
return NULL;
5663 for(
int jhit = 0; jhit < ntot; jhit++) {
5685 if(tmphitid == -1)
return NULL;
5686 refhit = cluster->
GetHit(tmphitid);
5706 double rc_of_min, rc_of_max;
5711 double u = chit->
GetU();
5712 double v = chit->
GetV();
5725 double rimin, rimax;
5726 r1 < r2 ? (rimin =
r1, rimax =
r2) : (rimin = r2, rimax = r1);
5728 rimin <
fRmin ? (rc_of_min = rc,
fRmin = rimin) : fRmin;
5729 rimax >
fRmax ? (rc_of_max = rc,
fRmax = rimax) : fRmax;
5738 double delta =
fabs(dv - du)/2.;
5739 du < dv ? (fUmin -= delta,
fUmax += delta) : (fVmin -= delta,
fVmax += delta);
5741 cout <<
"u_min " << fUmin <<
" u_max " << fUmax << endl;
5742 cout <<
"v_min " << fVmin <<
" v_max " <<
fVmax << endl;
5743 cout <<
"r_min " <<
fRmin <<
" r_max " <<
fRmax << endl;
5744 cout <<
"theta_min 0 theta_max 180" << endl;
5803 bool alreadythere =
false;
5808 if(
fVerbose > 1) cout <<
"\033[1;31m MAXPEAK " << maxpeak <<
", BREAK NOW! \033[0m" << endl;
5812 for(
size_t ialready = 0; ialready <
fFoundPeaks.size(); ialready++) {
5813 std::pair<double, double> foundthetar =
fFoundPeaks.at(ialready);
5814 double foundtheta = foundthetar.first;
5815 double foundr = foundthetar.second;
5817 if(theta_max == foundtheta && r_max == foundr) {
5820 alreadythere =
true;
5821 if(
fVerbose > 0) cout <<
"OH NO! THIS PEAK IS ALREADY THERE" << endl;
5826 if(alreadythere ==
false) {
5827 std::pair<double, double>
tr(theta_max, r_max);
5840 if(mode == 0 && alreadythere ==
true) {
5841 cout <<
"THIS PEAK IS ALREADY THERE" << endl;
5851 TMarker *mrk =
new TMarker(theta_max, r_max, 29);
5871 ycrot0 = 1 / (2 * fitp);
5872 xcrot0 = - fitm * ycrot0;
5873 R =
sqrt(xcrot0 * xcrot0 + ycrot0 * ycrot0);
5892 ycrot0 = 1/(2 * fita);
5893 xcrot0 = -fitb/(2 * fita);
5894 epsilon = -fitc * pow((1+(fitb*fitb)), -3/2);
5895 R = epsilon +
sqrt((xcrot0 * xcrot0)+(ycrot0 *ycrot0));
5926 fitp = 1 / (2 * ycrot0);
5927 fitm = - 2 * fitp * xcrot0;
5932 TObjArray sector[6];
5933 TObjArray *neighborings = NULL;
5954 neighborings =
new TObjArray();
5957 if(ihit == jhit)
continue;
5962 if(tube->
IsNeighboring(tubeID2) == kTRUE) neighborings->Add(hit2);
5996 for(
int jhit = 0; jhit < indiv.GetEntriesFast(); jhit++) {
6005 cout <<
" FILLHITMAP STARTING" << endl;
6020 cout <<
"CLUSTERIZATION <---------------" << endl;
6024 TObjArray neighborings;
6025 int clusterizedhits = 0;
6027 for(
int iseed = 0; iseed < seeds.GetEntriesFast(); iseed++) {
6032 if(seedhit->
IsUsed() == kTRUE)
continue;
6039 cluster->
AddHit(seedhit);
6045 cout <<
"SEED " << seedtubeID << endl;
6046 seedhit->
Draw(kRed);
6054 int nlastadded = 1, addedcounter = 0;
6057 while(nlastadded > 0) {
6076 for(
int iadded = nclusterhits - 1; iadded >= (nclusterhits - nlastadded); iadded--) {
6079 if(neighborings.GetEntriesFast() == 0)
continue;
6083 for(
int ineigh = 0; ineigh < neighborings.GetEntriesFast(); ineigh++)
6091 cluster->
AddHit(neighhit);
6099 neighhit->
Draw(kGreen);
6100 cout <<
"nof hits " << cluster->
GetNofHits() << endl;
6112 nlastadded = addedcounter;
6124 cout <<
"NOF TOTAL HITS " <<
stthitlist->
GetNofHits() <<
" NOF CLUSTERIZED HITS " << clusterizedhits << endl;
6130 cout <<
"we have " << candseeds.GetEntriesFast() <<
" canduidate seeds" << endl;
6131 for(
int jseed = 0; jseed < candseeds.GetEntriesFast(); jseed++) {
6141 cout <<
"CSEED " << endl;
6142 cseedhit->
Draw(kBlue);
6148 if(cseedhit->
IsUsed() == kTRUE) { cout <<
"already" << endl ;
continue; }
6150 int cseedtubeID = cseedhit->
GetTubeID();
6155 cluster->
AddHit(cseedhit);
6161 cout <<
"SEED " << cseedtubeID << endl;
6162 cseedhit->
Draw(kRed);
6170 int nlastadded = 1, addedcounter = 0;
6173 while(nlastadded > 0) {
6192 for(
int iadded = nclusterhits - 1; iadded >= (nclusterhits - nlastadded); iadded--) {
6195 if(neighborings.GetEntriesFast() == 0)
continue;
6199 for(
int ineigh = 0; ineigh < neighborings.GetEntriesFast(); ineigh++)
6207 cluster->
AddHit(neighhit);
6215 neighhit->
Draw(kGreen);
6216 cout <<
"nof hits " << cluster->
GetNofHits() << endl;
6228 nlastadded = addedcounter;
6268 int nofhitsinlay[30];
6269 for(
int ilay = 0; ilay < 30; ilay++) nofhitsinlay[ilay] = 0;
6287 int maxnoftracks = 1;
6292 for(
int ihit = 0; ihit < cluster.
GetNofHits(); ihit++) {
6299 if(nofhitsinlay[layid] <= 1)
continue;
6302 if(layid != tmplayid) {
6322 if(
counter1 == nofhitsinlay[layid])
continue;
6324 for(
int jhit = counter; jhit < counter + nofhitsinlay[layid] -
counter1; jhit++) {
6341 if(counter1 == nofhitsinlay[layid] - 1) {
6342 int noftracks = nofhitsinlay[layid] - isneigh;
6343 cout <<
"CLUSTER CONTAINS @ LAYER " << layid <<
" ACTUALLY " << nofhitsinlay[layid] <<
" - " << isneigh <<
" = " << noftracks <<
" TRACKS" << endl;
6344 if(noftracks > maxnoftracks) maxnoftracks = noftracks;
6350 cout <<
"THIS CLUSTER HAS A TOTAL OF " << maxnoftracks <<
" TRACKS" << endl;
6351 return maxnoftracks;
6387 cout <<
"COUNT TRACKS IN SKEW SECTOR" << endl;
6388 int nofhitsinlay[30];
6389 for(
int ilay = 0; ilay < 30; ilay++) nofhitsinlay[ilay] = 0;
6391 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++) {
6407 int maxnoftracks = 1;
6412 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++) {
6421 if(nofhitsinlay[layid] <= 1)
continue;
6424 if(layid != tmplayid) {
6444 if(
counter1 == nofhitsinlay[layid])
continue;
6446 for(
int jhit = counter; jhit < counter + nofhitsinlay[layid] -
counter1; jhit++) {
6463 if(counter1 == nofhitsinlay[layid] - 1) {
6464 int noftracks = nofhitsinlay[layid] - isneigh;
6465 cout <<
"CLUSTER CONTAINS @ LAYER " << layid <<
" ACTUALLY " << nofhitsinlay[layid] <<
" - " << isneigh <<
" = " << noftracks <<
" TRACKS" << endl;
6466 if(noftracks > maxnoftracks) maxnoftracks = noftracks;
6472 cout <<
"THIS CLUSTER HAS A TOTAL OF " << maxnoftracks <<
" TRACKS" << endl;
6473 return maxnoftracks;
6482 Double_t delta = 0, trasl[2] = {0., 0.};
6495 cout <<
"DELTA " << delta <<
" TRASL " << trasl[0] <<
" " << trasl[1] << endl;
6517 double theta_max, r_max;
6521 if(maxpeak < 4)
return NULL;
6527 TLine *line =
new TLine(-10.07, fitq + fitm * (-10.07), 10.07, fitq + fitm * (10.07));
6556 double rmin = R - R * 0.05;
6557 double rmax = R + R * 0.05;
6564 TArc *arcmin =
new TArc(xc, yc, rmin);
6565 TArc *arcmax =
new TArc(xc, yc, rmax);
6567 arcmin->SetFillStyle(0);
6568 arcmax->SetFillStyle(0);
6569 arcmin->SetLineColor(kGreen);
6570 arcmax->SetLineColor(kBlue);
6578 cout <<
"want to go to new cluster?" << endl;
6584 int startsecid = 1000, endsecid = -1, startlayid = 1000, endlayid = -1;
6585 double totaldistanceconf = 0;
6587 for(
int ihit = 0; ihit < cluster.
GetNofHits(); ihit++) {
6597 if(distance <= rmax && distance >= rmin) {
6598 thiscluster->
AddHit(hit);
6617 totaldistanceconf += distanceconf;
6629 thiscluster->
Draw(kRed);
6633 cout <<
"want to go to next cluster1?" << endl;
6638 if(startlayid != 0 || endlayid != 23) {
6645 if(distance <= rmax && distance >= rmin) {
6650 if(startsecid != 5 && endsecid != 5 && startsecid != 0 && endsecid != 0)
continue;
6670 cout <<
"want to go to next?" << endl;
6675 thiscluster->
AddHit(hit);
6686 thiscluster->
Draw(kRed);
6690 cout <<
"want to go to next cluster2?" << endl;
6707 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++)
6722 if(TMath::IsNaN(chit.
GetPosition().X()))
continue;
6728 mrk->SetMarkerColor(kRed);
6733 mrk2->SetMarkerColor(kRed);
6748 if(fitm == 0)
return kFALSE;
6756 cout <<
"wanna see the line?" << endl;
6757 TLine *line =
new TLine(-10.07, fitq + fitm * (-10.07), 10.07, fitq + fitm * (10.07));
6758 line->SetLineColor(2);
6783 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++)
6797 mrk->SetMarkerColor(kRed);
6805 mrk2->SetMarkerColor(kRed);
6819 cout <<
"previous " << xc <<
" " << yc <<
" " << R << endl;
6821 cout <<
"now " << xc <<
" " << yc <<
" " << R << endl;
6825 cout <<
"wanna see the line?" << endl;
6826 TLine *line =
new TLine(-10.07, fitp2 + fitm2 * (-10.07), 10.07, fitp2 + fitm2 * (10.07));
6827 line->SetLineColor(2);
6831 TArc *aline =
new TArc(xc, yc, R);
6832 aline->SetFillStyle(0);
6833 aline->SetLineColor(2);
6834 aline->Draw(
"SAME");
6859 for(
int ihit = 0; ihit < cluster->
GetNofHits(); ihit++)
6868 double sigma = 1e-5;
6870 if(hit->
IsGem()) sigma = 0.1;
6874 if(TMath::IsNaN(chit.
GetPosition().X()))
continue;
6880 mrk->SetMarkerColor(kRed);
6885 mrk2->SetMarkerColor(kRed);
6898 if(fita == 0)
return kFALSE;
6928 double xi = 0, yi = 0;
6930 fabs(yi1 - (fitm * xi1 + fitp)) <
fabs(yi2 - (fitm * xi2 + fitp)) ? (yi = yi1, xi = xi1) : (yi = yi2, xi = xi2);
6938 TVector2
vec(xc, yc);
6979 yb1 = vec.Y() +
sqrt(R * R - (xb1 - vec.X()) * (xb1 - vec.X()));
6982 yb2 = vec.Y() -
sqrt(R * R - (xb2 - vec.X()) * (xb2 - vec.X()));
6999 xb1 = (-(m*(q - vec.Y()) - vec.X()) +
sqrt((m*(q - vec.Y()) - vec.X())*(m*(q - vec.Y()) - vec.X()) - (m*m + 1)*((q - vec.Y())*(q - vec.Y()) + vec.X()*vec.X() - R *
R))) / (m*m + 1);
7002 xb2 = (-(m*(q - vec.Y()) - vec.X()) -
sqrt((m*(q - vec.Y()) - vec.X())*(m*(q - vec.Y()) - vec.X()) - (m*m + 1)*((q - vec.Y())*(q - vec.Y()) + vec.X()*vec.X() - R *
R))) / (m*m + 1);
7012 if(distb1 > distb2) xyb.Set(xb2, yb2);
7013 else xyb.Set(xb1, yb1);
7016 Double_t dist1 =
sqrt((xyb.Y() - y1)*(xyb.Y() - y1) + (xyb.X() - x1)*(xyb.X() - x1));
7017 Double_t dist2 =
sqrt((xyb.Y() - y2)*(xyb.Y() - y2) + (xyb.X() - x2)*(xyb.X() - x2));
7021 if(dist1 > dist2) xy =
new TVector2(x2, y2);
7022 else xy =
new TVector2(x1, y1);
7024 hit->
SetPosition(TVector3(xy->X(), xy->Y(), 0.0));
7041 if(hit->
IsSttSkew() == kFALSE)
continue;
7049 TVector3 first = tube->
GetPosition() + wireDirection * halflength;
7050 TVector3 second = tube->
GetPosition() - wireDirection * halflength;
7065 TVector2 intersection1, intersection2;
7068 if(nofintersections == 0)
continue;
7069 if(nofintersections >= 2) {
7070 cout <<
"ERROR: MORE THAN 1 INTERSECTION!!" << endl;
7077 TLine *l =
new TLine(first.X(), first.Y(), second.X(), second.Y());
7078 l->SetLineColor(kBlue);
7081 TMarker *mrk =
new TMarker(intersection1.X(), intersection1.Y(), 20);
7082 mrk->SetMarkerColor(kBlue);
7100 double beta = wireDirection.Phi();
7105 TVector2 rottangent(rtx, rty);
7106 rottangent = rottangent.Unit();
7117 Double_t rotm = rottangent.Y()/rottangent.X();
7127 x0a = (-rotp +
TMath::Sqrt(b * b + a * a * rotm * rotm)) / rotm;
7128 x0b = (-rotp -
TMath::Sqrt(b * b + a * a * rotm * rotm)) / rotm;
7131 double intxa = (x0a * b * b - rotm * rotp * a *
a) / (b * b + rotm * rotm * a * a);
7132 double intya = rotm * intxa + rotp;
7133 double intxb = (x0b * b * b - rotm * rotp * a *
a) / (b * b + rotm * rotm * a * a);
7134 double intyb = rotm * intxb + rotp;
7161 double y0a = y0anew;
7163 double y0b = y0bnew;
7169 TEllipse *ell1 =
new TEllipse(x0a, y0a, a, b, 0, 360, -beta);
7170 ell1->SetFillStyle(0);
7171 ell1->SetLineColor(4);
7173 TEllipse *ell2 =
new TEllipse(x0b, y0b, a, b, 0, 360, -beta);
7174 ell2->SetFillStyle(0);
7175 ell2->SetLineColor(6);
7178 TMarker *mrkinta =
new TMarker(intxa, intya, 20);
7179 mrkinta->SetMarkerColor(4);
7180 mrkinta->Draw(
"SAME");
7181 TMarker *mrkintb =
new TMarker(intxb, intyb, 20);
7182 mrkintb->SetMarkerColor(6);
7183 mrkintb->Draw(
"SAME");
7190 Double_t t = ((x0a + y0a) - (first.X() + first.Y())) / ((second.X() - first.X()) + (second.Y() - first.Y()));
7191 Double_t z0a = first.Z() + (second.Z() - first.Z()) * t;
7194 t = ((x0b + y0b) - (first.X() + first.Y())) / ((second.X() - first.X()) + (second.Y() - first.Y()));
7195 Double_t z0b = first.Z() + (second.Z() - first.Z()) * t;
7197 TVector3 center1(x0a, y0a, z0a);
7198 TVector3 center2(x0b, y0b, z0b);
7216 double dx = intxa - x0a;
7217 double dy = intya - y0a;
7218 TVector3 dxdy(dx, dy, 0.0);
7220 TVector3 tfirst = first + dxdy;
7221 TVector3 tsecond = second + dxdy;
7223 t = ((intxa + intya) - (tfirst.X() + tfirst.Y())) / ((tsecond.X() - tfirst.X()) + (tsecond.Y() - tfirst.Y()));
7224 double intza = tfirst.Z() + (tsecond.Z() - tfirst.Z()) * t;
7228 TLine *linezx1 =
new TLine(tfirst.X(), tfirst.Z(), tsecond.X(), tsecond.Z());
7229 linezx1->SetLineStyle(1);
7230 linezx1->Draw(
"SAME");
7231 TMarker *mrkza1 =
new TMarker(intxa, intza, 20);
7232 mrkza1->SetMarkerColor(kBlue - 9);
7233 mrkza1->Draw(
"SAME");
7237 tfirst = first - dxdy;
7238 tsecond = second - dxdy;
7240 t = ((intxb + intyb) - (tfirst.X() + tfirst.Y())) / ((tsecond.X() - tfirst.X()) + (tsecond.Y() - tfirst.Y()));
7241 double intzb = tfirst.Z() + (tsecond.Z() - tfirst.Z()) * t;
7243 TVector3 fin_intersection1(intxa, intya, intza);
7244 TVector3 fin_intersection2(intxb, intyb, intzb);
7263 double phi1 = track->
ComputePhi(fin_intersection1);
7264 double phi2 = track->
ComputePhi(fin_intersection2);
7266 PndTrkSkewHit *skewhit =
new PndTrkSkewHit(*hit, trackID, center1, fin_intersection1, phi1, center2, fin_intersection2, phi2, a, b, -1, beta);
7268 skewhitlist.
AddHit(skewhit);
7286 double phimin = 400, phimax = -1, zmin = 1000, zmax = -1;
7288 for(
int ihit = 0; ihit < skewhitlist->
GetNofHits(); ihit++) {
7290 if(!skewhit)
continue;
7294 double phi1 = skewhit->
GetPhi1();
7295 double phi2 = skewhit->
GetPhi2();
7297 if(phi1 < phimin) phimin = phi1;
7298 if(phi2 < phimin) phimin = phi2;
7299 if(fin_intersection1.Z() < zmin) zmin = fin_intersection1.Z();
7300 if(fin_intersection2.Z() < zmin) zmin = fin_intersection2.Z();
7302 if(phi1 > phimax) phimax = phi1;
7303 if(phi2 > phimax) phimax = phi2;
7304 if(fin_intersection1.Z() > zmax) zmax = fin_intersection1.Z();
7305 if(fin_intersection2.Z() > zmax) zmax = fin_intersection2.Z();
7320 hzphi->SetXTitle(
"#phi");
7321 hzphi->SetYTitle(
"#z");
7327 for(
int ihit = 0; ihit < skewhitlist->
GetNofHits(); ihit++) {
7329 if(!skewhit)
continue;
7334 double phi1 = skewhit->
GetPhi1();
7335 double phi2 = skewhit->
GetPhi2();
7337 TLine *linezphi =
new TLine(phi1, fin_intersection1.Z(), phi2, fin_intersection2.Z());
7339 linezphi->SetLineStyle(1);
7340 linezphi->Draw(
"SAME");
7342 TMarker *mrkzphi1 =
new TMarker(phi1, fin_intersection1.Z(), 20);
7345 mrkzphi1->SetMarkerColor(kBlue - 9);
7346 mrkzphi1->Draw(
"SAME");
7348 TMarker *mrkzphi2 =
new TMarker(phi2, fin_intersection2.Z(), 20);
7350 mrkzphi2->SetMarkerColor(kMagenta - 7);
7351 mrkzphi2->Draw(
"SAME");
7359 TH1F hz(
"hz",
"z", (zmax - zmin)/10., zmin, zmax);
7360 for(
int ihit = 0; ihit < skewhitlist->
GetNofHits(); ihit++) {
7362 if(!skewhit)
continue;
7365 hz.Fill(fin_intersection1.Z());
7366 hz.Fill(fin_intersection2.Z());
7368 int maxbinz = hz.GetMaximumBin();
7369 double mostprobZ = hz.GetBinCenter(maxbinz);
7375 for(
int ihit = 0; ihit < skewhitlist->
GetNofHits(); ihit++) {
7377 if(!skewhit)
continue;
7380 double phi1 = skewhit->
GetPhi1();
7381 double phi2 = skewhit->
GetPhi2();
7383 if(
fabs(fin_intersection1.Z() - mostprobZ) > 30. &&
fabs(fin_intersection2.Z() - mostprobZ) > 30.) {
7388 tmpskewhitlist.
AddHit(skewhit);
7390 tmpskewhitlist.
Sort();
7391 return tmpskewhitlist;
7397 if(
hzphi == NULL)
hzphi =
new TH2F(
"hzphi",
"z - phi plane", phimax - phimin, phimin, phimax, zmax - zmin, zmin, zmax);
7400 hzphi->GetXaxis()->SetLimits(phimin, phimax);
7401 hzphi->GetYaxis()->SetLimits(zmin, zmax);
7425 for (
size_t ihit = 0; ihit < pricand->
GetNHits(); ihit++) {
7430 std::vector< int > hitids = det_to_hitids[detId];
7431 if(std::find(hitids.begin(), hitids.end(), hitId) == hitids.end()) hitids.push_back(hitId);
7432 det_to_hitids[detId] = hitids;
7438 Int_t size = clref1.GetEntriesFast();
7443 size = clref2.GetEntriesFast();
7467 std::map< int, bool > hitidTousability;
7468 std::vector< int >
hits = det_to_hitids[detid];
7470 TClonesArray *array = NULL;
7476 for(
int ihit = 0; ihit < array->GetEntriesFast(); ihit++) {
7477 hitidTousability[ihit] =
true;
7478 if(find(hits.begin(), hits.end(), ihit) != hits.end()) hitidTousability[ihit] =
false;
7480 return hitidTousability;
void AddNeighboringsToHit(PndTrkHit *hit, TObjArray *hits)
PndTrkCluster CleanUpSkewHitList(PndTrkCluster *skewhitlist)
PndTrkClusterList CreateFullClusterization()
Bool_t SharedAt(PndTrkCluster *cluster2, double limit)
PndTrkCluster * fIndivCluster
PndTrkConformalTransform * conform
TClonesArray * fMvdStripHitArray
void AddTCA(Int_t detID, TClonesArray *array)
PndTrkSttHitList * stthitlist
PndTrkHit * GetHitByID(int id)
PndTrack ConvertToPndTrack()
TClonesArray * fTubeArray
TClonesArray * fPrimaryTrackArray
PndTrkHit * FindReferenceHit()
void Draw(Color_t color=kBlack)
friend F32vec4 sqrt(const F32vec4 &a)
int MergeTo(PndTrkCluster *cluster2)
PndTrkSttHitList * Instanciate()
void DeleteHitAndCompress(PndTrkHit *hit)
static T Sqrt(const T &x)
Double_t fStt_RealDistLimit
void AnalyticalFit2(PndTrkCluster *cluster, double fitm, double fitp, double &fitm2, double &fip2)
PndTrackCandHit GetSortedHit(UInt_t i)
void DrawNeighboringsToHit(PndTrkHit *hit)
PndTrkCluster GetCluster()
void Chi2Calculation(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t)
Double_t fMvdPix_RealDistLimit
PndTrkCluster * fIndivisibleHitList
Bool_t AnalyticalParabolaFit(PndTrkCluster *cluster, double xc, double yc, double R, double &fita, double &fitb, Double_t &fitc, Double_t &epsilon)
void SetTanL(double tanl)
PndTrkHit * FindMvdReferenceHit()
TClonesArray * fMvdPixelHitArray
PndTrkTrackList * fTrackList
TObjArray GetNeighboringsToHit(PndTrkHit *hit)
PndTrkGemHitList * gemhitlist
PndTrkConformalHitList * fConformalHitList
void SetPhi(Double_t phi)
int GetNofHitsInSector(int isec)
void SetSortVariable(Double_t sortvar)
PndTrkCluster * fFinalCluster
void AddNonCombiHits(Int_t detID, TClonesArray *array, std::map< int, bool > hitTousable)
void AddHit(PndTrkHit *hit)
void DrawConfHit(double x, double y, double r, int marker=2)
PndTrkHit * GetHitFromSector(int ihit, int isec)
Int_t ApplyLegendre(PndTrkCluster *cluster, double &theta_max, double &r_max)
void SetCluster(PndTrkCluster *cluster)
PndTrkCluster * CreateClusterAroundTrack(PndTrkTrack *track)
PndTrkHit * GetHit(int index)
PndTrkSkewHit * GetHit(int index)
TVector3 GetIntersection2()
PndTrkHit * GetHit(int index)
void SetOwnerValue(Bool_t enable=kTRUE)
void AddCluster(PndTrkCluster *cluster)
Double_t ComputePhi(TVector3 hit)
void RemoveHit(PndTrkHit *hit)
PndTrkSdsHitList * mvdpixhitlist
virtual InitStatus Init()
Bool_t MinuitFit(PndTrkCluster *cluster, double mstart, double qstart, double &fitm, double &fitq)
Bool_t MinuitFit2(PndTrkCluster *cluster, double xstart, double ystart, double rstart, double &xc, double &yc, double &R, double &sign)
std::map< int, bool > CombinatorialSuppression()
PndTrkCluster CreateSkewHitList(PndTrkTrack *track)
TClonesArray * fTrackCandArray
PndTrkGemHitList * Instanciate()
PndTrkSciTHitList * Instanciate()
void SetCenter(double x, double y)
Int_t RecreateHitArrays(std::map< int, std::vector< int > > &det_to_hitids)
void FromConformalToRealTrack(double fitm, double fitp, double &x0, double &y0, double &R)
void FillLegendreHisto(PndTrkCluster *cluster)
Int_t CountTracksInCluster(PndTrkCluster *cluster)
Double_t fit_distance(float x, float y, Double_t *par)
Bool_t IsNeighboring(int tubeID)
PndTrkTrack * GetTrack(Int_t index)
void DrawZGeometry(double phimin=0, double phimax=360, double zmin=-43, double zmax=113)
PndTrkSciTHitList * scithitlist
PndGeoSttPar * fSttParameters
void Clear(Option_t *="")
void AddTrack(PndTrkTrack *track)
PndTrackCand GetTrackCand()
void SetUsedFlag(Bool_t used)
PndTrkHit * FindMvdPixelReferenceHit()
void SetRadius(double radius)
PndTrkSdsHitList * InstanciateStrip()
virtual void Exec(Option_t *opt)
TClonesArray * fTrackArray
void FromRealToConformalTrack(double x0, double y0, double R, double &fitm, double &fitp)
static T ATan2(const T &y, const T &x)
void Draw(Color_t color=kBlack)
Int_t CountPossibleTracks()
void Apollonius(PndTrkCluster *cluster, std::vector< double > &X, std::vector< double > &Y, std::vector< double > &R)
void ComputeTraAndRot(PndTrkHit *hit, Double_t &delta, Double_t trasl[2])
Bool_t ParabolaFit(Double_t &fita, Double_t &fitb, Double_t &fitc)
void SetPhi2(Double_t phi2)
void SetPosition(TVector3 pos)
void DrawGeometryConf(double x1, double x2, double y1, double y2)
PndTrkLegendreTransform * legendre
FairTrackParP GetParamLast()
void CircleBy3Points(PndTrkHit *hit1, PndTrkHit *hit2, PndTrkHit *hit3, double &X, double &Y, double &R)
TClonesArray * fGemHitArray
TClonesArray * FillTubeArray()
void Chi2Calculation2(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t)
PndTrkNeighboringMap * fHitMap
std::map< int, bool > PrimaryCheck(Int_t detid, std::map< int, std::vector< int > > &det_to_hitids)
void ComputePlaneExtremities(PndTrkCluster *cluster)
char fMvdPixelBranch[200]
void Draw(Color_t color=kBlack)
PndTrkCluster * fSkewCluster
friend F32vec4 fabs(const F32vec4 &a)
PndTrkHit * FindSttReferenceHit(int isec=-1)
void AddHit(PndTrkSkewHit *shit)
void DrawTube(Color_t color)
static PndGeoHandling * Instance()
void DrawHits(PndTrkHitList *hitlist)
void RePrepareLegendre(PndTrkCluster *cluster)
Int_t ExtractLegendre(Int_t mode, double &theta_max, double &r_max)
Int_t ClusterToConformal(PndTrkCluster *cluster)
PndTrkCluster * fFinalSkewCluster
TClonesArray * fSciTHitArray
PndTrkSdsHitList * mvdstrhitlist
PndTrkTrack * LegendreFit(PndTrkCluster *cluster)
Bool_t DoesContain(PndTrkHit *hit)
TObjArray GetIndivisiblesToHit(PndTrkHit *hit)
TVector3 GetIntersection1()
std::vector< std::pair< double, double > > fFoundPeaks
Double_t GetXYDistance(PndTrkHit *fromhit)
Double_t fMvdPix_ConfDistLimit
TClonesArray * fSttPointArray
Bool_t StraightLineFit(Double_t &fitm, Double_t &fitp)
void IntersectionFinder(PndTrkHit *hit, double xc, double yc, double R)
void FromConformalToRealTrackParabola(double fita, double fitb, double fitc, double &x0, double &y0, double &R, double &epsilon)
PndSttMapCreator * fMapper
Double_t GetSortVariable()
void Clear(Option_t *opt="")
char fMvdStripBranch[200]
TClonesArray * fSttHitArray
Int_t FillConformalHitList(PndTrkCluster *cluster)
PndTrackCand * GetTrackCandPtr()
Double_t fMvdStr_RealDistLimit
TClonesArray * fTrkTrackArray
PndTrkGemCombinatorial * fCombiFinder
cout<<"the mcX, mcY, mcZ are "<< mcX<<","<< mcY<<","<< mcZ<< endl;vecmc.SetXYZ(mcX, mcY, mcZ);vecrc=hit-> GetPosition()
void SetPhi1(Double_t phi1)
Double_t fMvdStr_ConfDistLimit
PndTrkHit * FindMvdStripReferenceHit()
Bool_t AnalyticalFit(PndTrkCluster *cluster, double xc, double yc, double R, double &fitm, double &fitq)
TObjArray GetIndivisibles()
Int_t CountTracksInSkewSector(PndTrkCluster *cluster)
TVector3 GetWireDirection()
FairTrackParP GetParamFirst()
Double_t fStt_ConfDistLimit
void SetPointToFit(double x, double y, double sigma)
Double_t ComputePhiFrom(TVector3 hit, TVector3 from)
PndTrkSdsHitList * InstanciatePixel()