1 #include "TStopwatch.h"
4 #include "TClonesArray.h"
5 #include "basefunctions_ana.hxx"
24 gROOT->LoadMacro(
"$VMCWORKDIR/gconfig/rootlogon.C");
28 TString simFileName =
"Sim_Dpm_500.root";
29 TString parFileName =
"Sim_Dpm_500_params.root";
30 TString digiFileName =
"Sim_Dpm_500_digi.root";
31 TString recoFileName =
"Sim_Dpm_500_reco.root";
32 TString outFileName =
"Sim_Dpm_500_streamdisplay.root";
36 TFile filerecopixel(recoFileName.Data());
39 TClonesArray *sttsortedhits =
new TClonesArray(
"PndSttHit");
40 treedigi->SetBranchAddress(
"STTSortedHits",&sttsortedhits);
42 TTree *mvdstriprecotree = (TTree*)
filereco.Get(
"pndsim");
43 TClonesArray *mvdstripreco =
new TClonesArray(
"PndSdsHit");
44 mvdstriprecotree->SetBranchAddress(
"MVDHitsStrip",&mvdstripreco);
46 TTree *mvdpixelrecotree = (TTree*) filerecopixel.Get(
"pndsim");
47 TClonesArray *mvdpixelreco =
new TClonesArray(
"PndSdsHit");
48 mvdpixelrecotree->SetBranchAddress(
"MVDHitsPixel",&mvdpixelreco);
50 FairRunAna *
fRun=
new FairRunAna();
51 fRun->SetInputFile(simFileName.Data());
53 fRun->AddFriend(digiFileName.Data());
55 fRun->SetOutputFile(outFileName.Data());
57 FairRuntimeDb*
rtdb = fRun->GetRuntimeDb();
58 FairParRootFileIo*
parInput1 =
new FairParRootFileIo();
59 parInput1->open(parFileName.Data());
60 rtdb->setFirstInput(parInput1);
67 fRun->AddTask(sttFindTracks);
74 c1 =
new TCanvas(
"c1");
75 c1->Range(-42,-42,42,42);
76 c1->SetCanvasSize(1200, 1200);
90 Int_t timelimit = 2200;
92 Int_t timewindow = 245;
93 Int_t revtimewindow = 10;
94 Int_t mvdtimewindow = 20;
95 Int_t mvdrevtimewindow = 5;
99 Int_t stripnexttime = 0;
100 Int_t stripevpos = 0;
101 Int_t stripclonespos = 0;
102 Int_t pixelnexttime = 0;
103 Int_t pixelevpos = 0;
104 Int_t pixelclonespos = 0;
105 TObjArray *drawstackhits =
new TObjArray();
106 drawstackhits->SetName(
"drawstackhits");
107 drawstackhits->SetOwner(kTRUE);
108 TObjArray *drawstackstriphits =
new TObjArray();
109 drawstackstriphits->SetName(
"drawstackstriphits");
110 drawstackstriphits->SetOwner(kTRUE);
111 TObjArray *drawstackpixelhits =
new TObjArray();
112 drawstackpixelhits->SetName(
"drawstackpixelhits");
113 drawstackpixelhits->SetOwner(kTRUE);
116 TText* mytext =
new TText();
118 treedigi->GetEntry(0);
119 mvdstriprecotree->GetEntry(0);
120 mvdpixelrecotree->GetEntry(0);
125 TEllipse* mybglipse =
new TEllipse();
126 mybglipse->SetFillColor(kWhite);
127 Double_t xbg[8] = {2, 16, 16, 2, -2, -16, -16, -2};
128 Double_t ybg[8] = {17.2, 9.2, -9.2, -17.2, 17.2, 9.2, -9.2, -17.2};
129 TPolyLine* plinebg =
new TPolyLine(8, xbg, ybg);
130 plinebg->SetFillColor(kWhite);
134 Double_t yr[4] = {17.2, 9.2, -9.2, -17.2};
135 TPolyLine* pliner =
new TPolyLine(4, xr, yr);
137 pliner->SetLineColor(kBlack);
140 Double_t xl[4] = {-2, -16, -16, -2};
141 Double_t yl[4] = {17.2, 9.2, -9.2, -17.2};
142 TPolyLine* plinel =
new TPolyLine(4, xl, yl);
144 plinel->SetLineColor(kBlack);
147 Double_t xls[8] = {-2, -31.5, -31.5, -2, -2, -23.2, -23.2, -2};
148 Double_t yls[8] = {35.2, 18.5, -18.5, -35.2, -25.7, -13.5, 13.5, 25.7};
149 TPolyLine* plinels =
new TPolyLine(8, xls, yls);
151 Double_t xrs[8] = {2, 31.5, 31.5, 2, 2, 23.2, 23.2, 2};
152 Double_t yrs[8] = {35.2, 18.5, -18.5, -35.2, -25.7, -13.5, 13.5, 25.7};
153 TPolyLine* pliners =
new TPolyLine(8, xrs, yrs);
157 mybglipse->DrawEllipse(-2, 0, 39, 40.5, 90, 270, 0);
158 mybglipse->DrawEllipse(2, 0, 39, 40.5, 90, 270, 180);
164 mybglipse->DrawEllipse(0, 0, 14, 14, 0, 360, 0);
165 mybglipse->DrawEllipse(0, 0, 0.5, 0.5, 0, 360, 0);
167 for (Int_t currenttime = 0; currenttime < timelimit; currenttime+=5) {
169 GetHitsInTime(currenttime, timewindow, drawstackhits, nexttime, evpos, clonespos, treedigi, sttsortedhits);
170 GetHitsInTime(currenttime, mvdtimewindow, drawstackstriphits, stripnexttime, stripevpos, stripclonespos, mvdstriprecotree, mvdstripreco);
171 GetHitsInTime(currenttime, mvdtimewindow, drawstackpixelhits, pixelnexttime, pixelevpos, pixelclonespos, mvdpixelrecotree, mvdpixelreco);
174 cout << drawstackhits->GetEntriesFast() << endl;
175 if ((drawstackhits->IsEmpty()) && (drawstackstriphits->IsEmpty()) && (drawstackpixelhits->IsEmpty()))
continue;
178 RemoveOldHits(currenttime, revtimewindow, drawstackhits);
179 RemoveOldHits(currenttime, mvdrevtimewindow, drawstackstriphits);
180 RemoveOldHits(currenttime, mvdrevtimewindow, drawstackpixelhits);
185 mybglipse->DrawEllipse(-2, 0, 39, 40.5, 90, 270, 0);
186 mybglipse->DrawEllipse(2, 0, 39, 40.5, 90, 270, 180);
192 mybglipse->DrawEllipse(0, 0, 14, 14, 0, 360, 0);
193 mybglipse->DrawEllipse(0, 0, 0.5, 0.5, 0, 360, 0);
195 cout <<
" Drawing stuff " << endl;
197 DrawFairHits(drawstackstriphits);
198 DrawFairHits(drawstackpixelhits);
200 mytext->DrawText(-42, -42, TString::Format(
"%d ns",currenttime));
202 c1->Print(
"hitstream_tripletana.gif+15");
229 Double_t px = 0.5*(b*c*c-a*a*d-b*b*d+b*d*
d)/(b*c-a*d);
230 Double_t py = 0.5*(a*a*c+b*b*c-a*c*c-a*d*
d)/(b*c-a*d);
232 TEllipse mylipse(px, py, r, r);
233 mylipse.SetFillStyle(0);
234 mylipse.SetLineColor(linecolor);
239 Bool_t clockwise = (y1 > x1*y2/x2);
241 clockwise = !(clockwise);
243 cout <<
"Circle Parameters: " << px <<
" " << py <<
" " << x1 <<
" " << y1 <<
" " << x2 <<
" " << y2 <<
" " << TMath::RadToDeg()*zeroangle <<
" " << TMath::RadToDeg()*midangle <<
" " << TMath::RadToDeg()*endangle <<
" " << clockwise << endl;
257 mylipse.SetPhimin(TMath::RadToDeg()*endangle);
258 mylipse.SetPhimax(TMath::RadToDeg()*zeroangle);
259 mylipse.SetPhimin(180);
260 mylipse.SetPhimax(360);
274 mylipse.SetPhimin(TMath::RadToDeg()*zeroangle);
275 mylipse.SetPhimax(TMath::RadToDeg()*endangle);
276 mylipse.SetPhimin(0);
277 mylipse.SetPhimax(180);
280 mylipse.SetTheta(TMath::RadToDeg()*zeroangle);
284 void MakeTriplet(UInt_t leftID, TObjArray* stthits, TObjArray* stttubes) {
285 UInt_t neighbors = 0;
287 for (
int i = 0;
i < stthits->GetEntriesFast();
i++) {
288 rightID = ((
PndSttHit*)(stthits->At(
i)))->GetTubeID();
289 if (GetTubeDist(leftID, rightID, stttubes) < 1.5) {
296 UInt_t neighbors = 0;
298 for (
int i = 0;
i < stthits->GetEntriesFast();
i++) {
299 rightID = ((
PndSttHit*)(stthits->At(
i)))->GetTubeID();
301 if (GetTubeDist(leftID, rightID, stttubes) < 1.5) {
302 cout <<
"Neighbor found between tubes: " << leftID <<
", " << rightID <<
", " << GetTubeDist(leftID, rightID, stttubes) << endl;
306 return (neighbors > 2);
310 UInt_t neighbors = 0;
314 cmspoint->Set(cmsx, cmsy);
315 for (
int i = 0;
i < stthits->GetEntriesFast();
i++) {
316 rightID = ((
PndSttHit*)(stthits->At(
i)))->GetTubeID();
318 if (GetTubeDist(leftID, rightID, stttubes) < 1.5) {
319 cout <<
"Neighbor found between tubes: " << leftID <<
", " << rightID <<
", " << GetTubeDist(leftID, rightID, stttubes) << endl;
321 cmsx += ((
PndSttTube*)(stttubes->At(rightID)))->GetPosition()->X();
322 cmsy += ((
PndSttTube*)(stttubes->At(rightID)))->GetPosition()->Y();
325 cmspoint->Set(cmsx/neighbors, cmsy/neighbors);
326 return (neighbors > 2);
TFile filereco("MvdStt_Test_reco.root")
void MakeTriplet(UInt_t leftID, TObjArray *stthits, TObjArray *stttubes)
static T Sqrt(const T &x)
PndSttTrackFinderReal * sttTrackFinder
TEllipse DrawOriginTrackLipse(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Color_t linecolor=kBlack)
void DrawIsochrones(Int_t currenttime, TObjArray *stthits, TObjArray *stttubes, Int_t *timeviolationcount)
Bool_t HasNeighbor(UInt_t leftID, TObjArray *stthits, TObjArray *stttubes)
static T ATan2(const T &y, const T &x)
TFile filedigi("testdigi.root")
PndSttFindTracks * sttFindTracks
Double_t ATan2Pi(Double_t y, Double_t x)
TClonesArray * FillTubeArray()
Bool_t HasNeighborCMS(UInt_t leftID, TObjArray *stthits, TObjArray *stttubes, TVector2 *cmspoint)
FairParRootFileIo * parInput1
void AddHitCollectionName(char *hitCollectionName, char *pointCollectionName)