FairRoot/PandaRoot
Functions
runTripletAna.C File Reference
#include "TStopwatch.h"
#include "TTree.h"
#include "TFile.h"
#include "TClonesArray.h"
#include "basefunctions_ana.hxx"

Go to the source code of this file.

Functions

int runTripletAna ()
 
TEllipse DrawOriginTrackLipse (Double_t x1, Double_t y1, Double_t x2, Double_t y2, Color_t linecolor=kBlack)
 
void MakeTriplet (UInt_t leftID, TObjArray *stthits, TObjArray *stttubes)
 
Bool_t HasNeighbor (UInt_t leftID, TObjArray *stthits, TObjArray *stttubes)
 
Bool_t HasNeighborCMS (UInt_t leftID, TObjArray *stthits, TObjArray *stttubes, TVector2 *cmspoint)
 
Double_t ATan2Pi (Double_t y, Double_t x)
 

Function Documentation

Double_t ATan2Pi ( Double_t  y,
Double_t  x 
)

Definition at line 331 of file runTripletAna.C.

References angle, CAMath::ATan2(), Double_t, and Pi.

Referenced by DrawOriginTrackLipse().

331  {
332  Double_t angle = TMath::ATan2(y, x);
333  if (angle < 0) {
334  angle = 2*TMath::Pi() + angle;
335  }
336  return angle;
337 }
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
CAMath::ATan2
static T ATan2(const T &y, const T &x)
x
Double_t x
Definition: createRootFscGeometryFile.C:159
y
Double_t y
Definition: createRootFscGeometryFile.C:159
angle
Double_t angle
Definition: createROOTEnDrcGeo.C:13
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
TEllipse DrawOriginTrackLipse ( Double_t  x1,
Double_t  y1,
Double_t  x2,
Double_t  y2,
Color_t  linecolor = kBlack 
)

Definition at line 224 of file runTripletAna.C.

References a, ATan2Pi(), b, Bool_t, c, d, Double_t, Pi, r, and CAMath::Sqrt().

224  {
225  Double_t a= x1;
226  Double_t b= y1;
227  Double_t c= x2;
228  Double_t d= y2;
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);
231  Double_t r = TMath::Sqrt(px*px+py*py);
232  TEllipse mylipse(px, py, r, r);
233  mylipse.SetFillStyle(0);
234  mylipse.SetLineColor(linecolor);
235  Double_t zeroangle = ATan2Pi(-py, -px);
236  Double_t midangle = ATan2Pi(y1-py, x1-px);
237  Double_t endangle = ATan2Pi(y2-py, x2-px);
238  Double_t mirrorangle = zeroangle - TMath::Pi();
239  Bool_t clockwise = (y1 > x1*y2/x2);
240  if (x2 < 0) {
241  clockwise = !(clockwise);
242  }
243  cout << "Circle Parameters: " << px << " " << py << " " << x1 << " " << y1 << " " << x2 << " " << y2 << " " << TMath::RadToDeg()*zeroangle << " " << TMath::RadToDeg()*midangle << " " << TMath::RadToDeg()*endangle << " " << clockwise << endl;
244  if (clockwise) {
245  /*
246  Double_t diffangle = zeroangle - endangle;
247  if (diffangle < 0) {
248  diffangle = zeroangle + 2*TMath::Pi() - endangle;
249  }
250  if (diffangle < TMath::Pi()) {
251  endangle = zeroangle - TMath::Pi();
252  if (endangle < 0) {
253  endangle += 2*TMath::Pi();
254  }
255  }
256  */
257  mylipse.SetPhimin(TMath::RadToDeg()*endangle);
258  mylipse.SetPhimax(TMath::RadToDeg()*zeroangle);
259  mylipse.SetPhimin(180);
260  mylipse.SetPhimax(360);
261  } else {
262  /*
263  Double_t diffangle = endangle - zeroangle;
264  if (diffangle < 0) {
265  diffangle = 2*TMath::Pi() - zeroangle + endangle;
266  }
267  if (diffangle < TMath::Pi()) {
268  endangle = zeroangle + TMath::Pi();
269  if (endangle > 2*TMath::Pi()) {
270  endangle -= 2*TMath::Pi();
271  }
272  }
273  */
274  mylipse.SetPhimin(TMath::RadToDeg()*zeroangle);
275  mylipse.SetPhimax(TMath::RadToDeg()*endangle);
276  mylipse.SetPhimin(0);
277  mylipse.SetPhimax(180);
278  }
279 
280  mylipse.SetTheta(TMath::RadToDeg()*zeroangle);
281  return mylipse;
282 }
r
double r
Definition: RiemannTest.C:14
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
a
Int_t a
Definition: anaLmdDigi.C:126
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
ATan2Pi
Double_t ATan2Pi(Double_t y, Double_t x)
Definition: runTripletAna.C:331
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
Bool_t HasNeighbor ( UInt_t  leftID,
TObjArray *  stthits,
TObjArray *  stttubes 
)

Definition at line 295 of file runTripletAna.C.

References i.

295  {
296  UInt_t neighbors = 0;
297  UInt_t rightID = 0;
298  for (int i = 0; i < stthits->GetEntriesFast(); i++) {
299  rightID = ((PndSttHit*)(stthits->At(i)))->GetTubeID();
300  //cout << "Checking distance between tubes: " << leftID << ", " << rightID << endl;
301  if (GetTubeDist(leftID, rightID, stttubes) < 1.5) {
302  cout << "Neighbor found between tubes: " << leftID << ", " << rightID << ", " << GetTubeDist(leftID, rightID, stttubes) << endl;
303  neighbors++;
304  }
305  }
306  return (neighbors > 2);
307 }
i
Int_t i
Definition: run_full.C:25
Bool_t HasNeighborCMS ( UInt_t  leftID,
TObjArray *  stthits,
TObjArray *  stttubes,
TVector2 *  cmspoint 
)

Definition at line 309 of file runTripletAna.C.

References Double_t, and i.

309  {
310  UInt_t neighbors = 0;
311  UInt_t rightID = 0;
312  Double_t cmsx = 0;
313  Double_t cmsy = 0;
314  cmspoint->Set(cmsx, cmsy);
315  for (int i = 0; i < stthits->GetEntriesFast(); i++) {
316  rightID = ((PndSttHit*)(stthits->At(i)))->GetTubeID();
317  //cout << "Checking distance between tubes: " << leftID << ", " << rightID << endl;
318  if (GetTubeDist(leftID, rightID, stttubes) < 1.5) {
319  cout << "Neighbor found between tubes: " << leftID << ", " << rightID << ", " << GetTubeDist(leftID, rightID, stttubes) << endl;
320  neighbors++;
321  cmsx += ((PndSttTube*)(stttubes->At(rightID)))->GetPosition()->X();
322  cmsy += ((PndSttTube*)(stttubes->At(rightID)))->GetPosition()->Y();
323  }
324  }
325  cmspoint->Set(cmsx/neighbors, cmsy/neighbors);
326  return (neighbors > 2);
327 }
i
Int_t i
Definition: run_full.C:25
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
void MakeTriplet ( UInt_t  leftID,
TObjArray *  stthits,
TObjArray *  stttubes 
)

Definition at line 284 of file runTripletAna.C.

References i.

284  {
285  UInt_t neighbors = 0;
286  UInt_t rightID = 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) {
290  neighbors++;
291  }
292  }
293 }
i
Int_t i
Definition: run_full.C:25
int runTripletAna ( )

Definition at line 17 of file runTripletAna.C.

References PndSttFindTracks::AddHitCollectionName(), c1, Double_t, DrawIsochrones(), filedigi(), filereco(), PndSttMapCreator::FillTubeArray(), fRun, iVerbose, parInput1, rootlogon(), rtdb, sttFindTracks, sttTrackFinder, timer, treedigi, and TString.

17  {
18  gROOT->Reset();
19  // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
20  Int_t iVerbose = 0;
21 
22  TStopwatch timer;
23  timer.Start();
24  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
25  rootlogon();
26  //gSystem->Load("libSttMvdTracking");
27 
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";
33 
34  TFile filedigi(digiFileName.Data());
35  TFile filereco(recoFileName.Data());
36  TFile filerecopixel(recoFileName.Data());
37 
38  TTree *treedigi = (TTree*) filedigi.Get("pndsim");
39  TClonesArray *sttsortedhits = new TClonesArray("PndSttHit");
40  treedigi->SetBranchAddress("STTSortedHits",&sttsortedhits);
41 
42  TTree *mvdstriprecotree = (TTree*) filereco.Get("pndsim");
43  TClonesArray *mvdstripreco = new TClonesArray("PndSdsHit");
44  mvdstriprecotree->SetBranchAddress("MVDHitsStrip",&mvdstripreco);
45 
46  TTree *mvdpixelrecotree = (TTree*) filerecopixel.Get("pndsim");
47  TClonesArray *mvdpixelreco = new TClonesArray("PndSdsHit");
48  mvdpixelrecotree->SetBranchAddress("MVDHitsPixel",&mvdpixelreco);
49 
50  FairRunAna *fRun= new FairRunAna();
51  fRun->SetInputFile(simFileName.Data());
52  //fRun->AddFriend(recoFile.Data());
53  fRun->AddFriend(digiFileName.Data());
54  //fRun->AddFriend(trackF.Data());
55  fRun->SetOutputFile(outFileName.Data());
56  //fRun->RunWithTimeStamps();
57  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
58  FairParRootFileIo* parInput1 = new FairParRootFileIo();
59  parInput1->open(parFileName.Data());
60  rtdb->setFirstInput(parInput1);
61 
62  //works around the "geometry not supported by map" error
63  PndSttTrackFinderReal* sttTrackFinder = new PndSttTrackFinderReal(0);
64  PndSttFindTracks* sttFindTracks = new PndSttFindTracks("Track Finder", "FairTask", sttTrackFinder, iVerbose);
65  sttFindTracks->AddHitCollectionName("STTHit", "STTPoint");
66  //sttFindTracks->SetPersistence(kFALSE);
67  fRun->AddTask(sttFindTracks);
68  fRun->Init();
69 
70  PndGeoSttPar* fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar");
71  PndSttMapCreator* mapper = new PndSttMapCreator(fSttParameters);
72  TClonesArray* tubearray = mapper->FillTubeArray();
73 
74  c1 = new TCanvas("c1");
75  c1->Range(-42,-42,42,42);
76  c1->SetCanvasSize(1200, 1200);
77  //TEllipse mylipse(0, 0, 0.5, 0.5);
78 
79  // single straw tube simulation -----------------------
80  //PndSttSingleStraw stt;
81  //setting the single straw tube simulation constants
82  // 3 options currently available:
83  // TConst(tube radius (cm), gas pressure (bar), Ar%, CO2%)
84  // stt.TConst(0.4, 1, 0.9, 0.1);
85  // stt.TConst(0.5, 1, 0.9, 0.1);
86  //stt.TConst(0.5, 2, 0.8, 0.2);
87  // wire positioning
88  //stt.PutWireXYZ(0., 0., -75., 0., 0., 75.);
89 
90  Int_t timelimit = 2200;
91  //Int_t timelimit = 30;
92  Int_t timewindow = 245;
93  Int_t revtimewindow = 10;
94  Int_t mvdtimewindow = 20;
95  Int_t mvdrevtimewindow = 5;
96  Int_t nexttime = 0;
97  Int_t evpos = 0;
98  Int_t clonespos = 0;
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);
114  //PndSttHit* stthit = 0;
115  //PndSttTube* tube = 0;
116  TText* mytext = new TText();
117 
118  treedigi->GetEntry(0);
119  mvdstriprecotree->GetEntry(0);
120  mvdpixelrecotree->GetEntry(0);
121 
122  //Overlapping background goes to loops when canvas is cleared
123  //mylipse.DrawEllipse(-2, 0, 39, 40.5, 90, 270, 0);
124  //mylipse.DrawEllipse(2, 0, 39, 40.5, 90, 270, 180);
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);
131  //plinebg->SetLineColor(kWhite);
132  //plinebg->Draw("f");
133  Double_t xr[4] = {2, 16, 16, 2};
134  Double_t yr[4] = {17.2, 9.2, -9.2, -17.2};
135  TPolyLine* pliner = new TPolyLine(4, xr, yr);
136  //pliner->SetFillColor(kWhite);
137  pliner->SetLineColor(kBlack);
138  //pliner->Draw("f");
139  //pliner->Draw();
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);
143  //plinel->SetFillColor(kWhite);
144  plinel->SetLineColor(kBlack);
145  //plinel->Draw("f");
146  //plinel->Draw();
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);
150  //plinels->Draw();
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);
154  //pliners->Draw();
155 
156  //Draw Background
157  mybglipse->DrawEllipse(-2, 0, 39, 40.5, 90, 270, 0);
158  mybglipse->DrawEllipse(2, 0, 39, 40.5, 90, 270, 180);
159  plinebg->Draw("f");
160  pliner->Draw();
161  plinel->Draw();
162  pliners->Draw();
163  plinels->Draw();
164  mybglipse->DrawEllipse(0, 0, 14, 14, 0, 360, 0);
165  mybglipse->DrawEllipse(0, 0, 0.5, 0.5, 0, 360, 0);
166 
167  for (Int_t currenttime = 0; currenttime < timelimit; currenttime+=5) {
168  //Get events up to mytime + timewindow
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);
172 
173 
174  cout << drawstackhits->GetEntriesFast() << endl;
175  if ((drawstackhits->IsEmpty()) && (drawstackstriphits->IsEmpty()) && (drawstackpixelhits->IsEmpty())) continue;
176 
177  //Discard events < mytime
178  RemoveOldHits(currenttime, revtimewindow, drawstackhits);
179  RemoveOldHits(currenttime, mvdrevtimewindow, drawstackstriphits);
180  RemoveOldHits(currenttime, mvdrevtimewindow, drawstackpixelhits);
181 
182  //Calculate Isochrones and Draw
183  c1->Clear();
184  c1->cd();
185  mybglipse->DrawEllipse(-2, 0, 39, 40.5, 90, 270, 0);
186  mybglipse->DrawEllipse(2, 0, 39, 40.5, 90, 270, 180);
187  plinebg->Draw("f");
188  pliner->Draw();
189  plinel->Draw();
190  pliners->Draw();
191  plinels->Draw();
192  mybglipse->DrawEllipse(0, 0, 14, 14, 0, 360, 0);
193  mybglipse->DrawEllipse(0, 0, 0.5, 0.5, 0, 360, 0);
194 
195  cout << " Drawing stuff " << endl;
196  DrawIsochrones(currenttime, drawstackhits, tubearray);
197  DrawFairHits(drawstackstriphits);
198  DrawFairHits(drawstackpixelhits);
199 
200  mytext->DrawText(-42, -42, TString::Format("%d ns",currenttime));
201  c1->Update();
202  c1->Print("hitstream_tripletana.gif+15");
203  }
204  //c1->Print("hitstream.gif++");
205 
206 /*
207  for (int i = 0; i < tubearray->GetEntriesFast(); i++) {
208  tube = (PndSttTube*)(tubearray->At(i));
209  if (!tube) continue;
210  mybglipse.SetX1(tube->GetPosition().X());
211  mybglipse.SetY1(tube->GetPosition().Y());
212  mybglipse.SetR1(0.5);
213  mybglipse.SetR2(0.5);
214  if ( ((i > 104) && (i < 215))
215  || ((i > 714) && (i < 855))
216  || ((i > 2956) && (i < 3175)) )
217  mybglipse.DrawClone();
218  std::cout << i << " " << tube->GetPosition().X() << " " << tube->GetPosition().Y() << " " << tube->GetPosition().Z() << " " << tube->GetWireDirection().X() << " " << tube->GetWireDirection().Y() << " " << tube->GetWireDirection().Z() << std::endl;
219  }
220 */
221  return 0;
222 }
filereco
TFile filereco("MvdStt_Test_reco.root")
sttTrackFinder
PndSttTrackFinderReal * sttTrackFinder
Definition: run_tracking_sttcombi.C:64
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
DrawIsochrones
void DrawIsochrones(Int_t currenttime, TObjArray *stthits, TObjArray *stttubes, Int_t *timeviolationcount)
Definition: runTimestampCheck.C:163
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
treedigi
TTree * treedigi
Definition: checkdedx_helixhit.C:19
timer
TStopwatch timer
Definition: hit_dirc.C:51
filedigi
TFile filedigi("testdigi.root")
sttFindTracks
PndSttFindTracks * sttFindTracks
Definition: run_tracking_sttcombi.C:65
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
PndSttMapCreator::FillTubeArray
TClonesArray * FillTubeArray()
Definition: PndSttMapCreator.cxx:128
rootlogon
rootlogon()
Definition: outdated/mpiTools/macros/emc/rootlogon.C:1
c1
c1
Definition: plot_dirc.C:35
parInput1
FairParRootFileIo * parInput1
Definition: hit_dirc.C:67
iVerbose
Int_t iVerbose
Definition: Lars/runMvdDigi.C:7
PndSttFindTracks::AddHitCollectionName
void AddHitCollectionName(char *hitCollectionName, char *pointCollectionName)
Definition: PndSttFindTracks.cxx:151
Generated on Wed Apr 3 2019 08:44:23 for FairRoot/PandaRoot by   doxygen 1.8.5