1 #include "drawhelpers.hxx"
5 if (storagefile.Get(currenthisto->GetName()) != 0) {
6 currenthisto->Add((TH1*)storagefile.Get(currenthisto->GetName()));
8 cout << currenthisto->GetName() <<
": " << currenthisto->GetEntries() << endl;
9 currenthisto->Write(currenthisto->GetName(), TObject::kWriteDelete);
18 kDCH,
kDRC,
kDSK,
kEMC,
kGEM,
kLUMI,
kMDT,
kMVD,
kRPC,
kSTT,
kTPC,
kTOF,
kFTS,
kHYPG,
kHYP};
22 gROOT->LoadMacro(
"$VMCWORKDIR/gconfig/rootlogon.C");
25 TDatabasePDG *pdg =
new TDatabasePDG();
30 TString simFileName =
"Sim_Dpm_500.root";
31 TString parFileName =
"Sim_Dpm_500_params.root";
32 TString digiFileName =
"Sim_Dpm_500_digi.root";
33 TString recoFileName =
"Sim_Dpm_500_reco.root";
34 TString outFileName =
"Sim_Dpm_500_streamdisplay.root";
35 TString histofilename =
"onlinetrkmccheck.root";
39 TFile filerecopixel(recoFileName.Data());
42 TTree *recotree = (TTree*)
filereco.Get(
"pndsim");
44 TFile simfile(simFileName.Data());
45 TTree* mctree = (TTree*) simfile.Get(
"pndsim");
46 TClonesArray* mcarray =
new TClonesArray(
"PndMCTrack");
47 mctree->SetBranchAddress(
"MCTrack", &mcarray);
49 FairRunAna *
fRun =
new FairRunAna();
50 fRun->SetInputFile(simFileName.Data());
51 fRun->AddFriend(recoFileName.Data());
52 fRun->AddFriend(digiFileName.Data());
53 fRun->SetOutputFile(outFileName.Data());
55 FairRuntimeDb*
rtdb = fRun->GetRuntimeDb();
57 PndOnlineGeometryManager *online_geometry =
new PndOnlineGeometryManager(rtdb, parFileName.Data());
59 PndOnlineManager *online =
new PndOnlineManager();
60 online->AddGeometryManager(online_geometry);
63 online->LoadStream(PndOnline::kMVDPixel,
"MVDHitsPixel",stthitlifetime);
64 online->LoadStream(PndOnline::kMVDStrip,
"MVDHitsStrip",stthitlifetime);
67 TClonesArray* tubearray = online_geometry->GetDetectorGeometry(
PndOnline::kSTT);
70 PndOnlineSttTripletFinder *triplet_finder =
new PndOnlineSttTripletFinder(online, tubearray);
71 online->AddTask((FairTask*)triplet_finder);
74 TCanvas*
c1 =
new TCanvas(
"c1");
75 c1->Range(-42,-42,42,42);
76 c1->SetCanvasSize(1200, 1200);
77 TText* mytext =
new TText();
84 int final_time = maximumTime;
86 online->SetHESRRevolutionDuration(delta_t);
90 TObjArray* online_fairhits_stt = 0;
91 TObjArray* online_tracks = 0;
92 TObjArray* online_tracks2 = 0;
93 TObjArray* online_fairhits_mvd = 0;
95 std::vector<double> mctracktime;
96 std::vector<double> mctrackpt;
97 std::vector<double> mctrackphi;
98 std::vector<int> mctrackstatus;
99 std::vector<double> recotracktime;
100 std::vector<double> recotrackpt;
101 std::vector<double> recotrackphi;
102 std::vector<int> recotrackstatus;
103 std::vector<int> recotrackmcid;
104 int totalmctracks = 0;
105 int goodmctracks = 0;
106 int totalrecotracks = 0;
107 int goodrecotracks = 0;
108 int norecotracks = 0;
109 int deltaphifailed = 0;
110 int deltaptfailed = 0;
111 int deltabothfailed = 0;
113 int globalmaxmcid = 0;
123 TH2D* hmcall =
new TH2D(
"hmcall",
"hmcall", 1500, -15, 15, 75, 0, 1.5);
124 TH2D* hmcprime =
new TH2D(
"hmcprime",
"hmcprime", 1500, -15, 15, 75, 0, 1.5);
125 TH2D* hmcrecoable =
new TH2D(
"hmcrecoable",
"hmcrecoable", 1500, -15, 15, 75, 0, 1.5);
126 TH2D* hdafuq =
new TH2D(
"hdafuq",
"hdafuq", 1500, -15, 15, 75, 0, 1.5);
127 TH2D* hmcrecoed =
new TH2D(
"hmcrecoed",
"hmcrecoed", 1500, -15, 15, 75, 0, 1.5);
128 TH2D* hmcnotrecoed =
new TH2D(
"hmcnotrecoed",
"hmcnotrecoed", 1500, -15, 15, 75, 0, 1.5);
130 TH1D* hphires =
new TH1D(
"hphires",
"hphires", 3000, -15, 15);
131 TH2D* hptres =
new TH2D(
"hptres",
"hptres", 300, -1500, 1500, 200, -1, 1);
132 TH1D* hphiresmc =
new TH1D(
"hphiresmc",
"hphiresmc", 3000, -15, 15);
133 TH2D* hptresmc =
new TH2D(
"hptresmc",
"hptresmc", 300, -1500, 1500, 200, -1, 1);
134 TH1D* hphiresmc90 =
new TH1D(
"hphiresmc90",
"hphiresmc90", 3000, -15, 15);
135 TH2D* hptresmc90 =
new TH2D(
"hptresmc90",
"hptresmc90", 300, -1500, 1500, 200, -1, 1);
136 TH1D* hphiresmc80 =
new TH1D(
"hphiresmc80",
"hphiresmc80", 3000, -15, 15);
137 TH2D* hptresmc80 =
new TH2D(
"hptresmc80",
"hptresmc80", 300, -1500, 1500, 200, -1, 1);
138 TH1D* hphiresmcbad =
new TH1D(
"hphiresmcbad",
"hphiresmcbad", 3000, -15, 15);
139 TH2D* hptresmcbad =
new TH2D(
"hptresmcbad",
"hptresmcbad", 300, -1500, 1500, 200, -1, 1);
142 for (current_time = 0; current_time < final_time; current_time += delta_t) {
145 online->ClearTracks();
146 online->LoadHits( delta_t );
151 online_tracks = online->GetTrackObjectList();
155 if (online_tracks->GetEntriesFast() > 0) {
156 PndOnlineTrack* firstTrack = (PndOnlineTrack*) online_tracks->At(firsti);
157 while (firstTrack->DrawMode() != PndOnlineTrack::kVertexMomentumAsCircleHack) {
159 if (firsti > online_tracks->GetEntriesFast()) {
160 cout <<
"No drawable online Track in Array. We have a severe problem." << endl;
163 firstTrack = (PndOnlineTrack*) online_tracks->At(firsti);
165 maxmcid = firstTrack->MCEntryNumber();
166 minmcid = firstTrack->MCEntryNumber();
168 for (
int i = 0;
i < online_tracks->GetEntriesFast();
i++) {
169 PndOnlineTrack* myTrack = (PndOnlineTrack*) online_tracks->At(
i);
170 if (myTrack->DrawMode() == PndOnlineTrack::kVertexMomentumAsCircleHack) {
172 innerhit.Set(myTrack->Vertex().x(), myTrack->Vertex().y());
174 outerhit.Set(myTrack->Momentum().x(), myTrack->Momentum().y());
175 TVector2 circleorigin;
178 CalculateCircle(&innerhit, &outerhit, &circleorigin, &circleradius, &clockwise);
179 Double_t trackpt = circleradius * 2.0 * 3.0;
180 Double_t orthophi = circleorigin.Phi();
189 cout <<
"Reco Track MCEntryNumber: " << myTrack->MCEntryNumber() << endl;
190 cout <<
"Reco Track MCEntryMajority: " << myTrack->MCEntryMajority() << endl;
191 cout <<
"Reco Track MCTrackNumber: " << myTrack->MCTrackNumber() << endl;
192 cout <<
"Reco Track MCTrackMajority: " << myTrack->MCTrackMajority() << endl;
193 cout <<
"Reco Track HitCount: " << myTrack->HitCount() << endl;
195 mctree->GetEntry(myTrack->MCEntryNumber());
196 cout <<
"MC Entry count: " << mcarray->GetEntriesFast() << endl;
206 if (calcphimc < 0) calcphimc += 2*
TMath::Pi();
208 if (calcphireco < 0) calcphireco += 2*
TMath::Pi();
209 Double_t deltaphi = calcphireco - calcphimc;
217 Int_t chargesign = 1;
218 if (pdg->GetParticle(mctrack->
GetPdgCode())->Charge() < 0) {
221 hphiresmc->Fill(deltaphi*chargesign);
222 hptresmc->Fill(mcpt, deltapt/mcpt);
223 if ( (myTrack->MCEntryMajority() > 0.9) && (myTrack->MCTrackMajority() > 0.9) ) {
224 hphiresmc90->Fill(deltaphi*chargesign);
225 hptresmc90->Fill(mcpt, deltapt/mcpt);
227 if ( (myTrack->MCEntryMajority() > 0.8) && (myTrack->MCTrackMajority() > 0.8) ) {
228 hphiresmc80->Fill(deltaphi*chargesign);
229 hptresmc80->Fill(mcpt, deltapt/mcpt);
231 hphiresmcbad->Fill(deltaphi*chargesign);
232 hptresmcbad->Fill(mcpt, deltapt/mcpt);
236 cout <<
"MC Track has NULL pointer. Oooops." << endl;
241 if (mctrack->
GetMomentum().Pt() < 0.052) isrecoable =
false;
243 cout <<
"Low pt hit" << endl;
245 if ( (myTrack->MCEntryMajority() > 0.9) && (myTrack->MCTrackMajority() > 0.9) ) {
249 if ( (myTrack->MCEntryMajority() > 0.8) && (myTrack->MCTrackMajority() > 0.8) ) {
256 int mcentrynumber = myTrack->MCEntryNumber();
257 if (mcentrynumber > maxmcid) {
258 maxmcid = mcentrynumber;
260 if (mcentrynumber < minmcid) {
261 minmcid = mcentrynumber;
269 recotrackpt.push_back(trackpt);
270 recotrackphi.push_back(trackphi);
271 recotrackstatus.push_back(0);
272 recotrackmcid.push_back(mcentrynumber);
276 if (maxmcid > globalmaxmcid) {
277 globalmaxmcid = maxmcid;
280 for (
int j = minmcid; j < mctree->GetEntriesFast(); j++) {
281 cout <<
"MC outer loop: " << j << endl;
282 if (j > globalmaxmcid) {
286 cout <<
"MC Entry count: " << mcarray->GetEntriesFast() << endl;
287 for (
int k = 0; k < mcarray->GetEntriesFast(); k++) {
290 cout <<
"MC Track has NULL pointer. Oooops." << endl;
294 bool isrecoable =
false;
299 cout <<
"Critical Track at Event: " << j <<
" Track " << k <<
" Particle: " << mctrack->
GetPdgCode() << endl;
310 if (mctrack->
GetMomentum().Pt() < 0.052) isrecoable =
false;
312 cout <<
"Low pt hit" << endl;
318 if (pdg->GetParticle(mctrack->
GetPdgCode())->
Mass() > 1) isrecoable =
false;
321 cout <<
"Low betagamma particle" << endl;
325 cout <<
"Meh, reco possibility mismatch." << endl;
327 cout <<
"Dafuq at Event: " << j <<
" Track " << k <<
" Particle: " << mctrack->
GetPdgCode() <<
" Charge: " << pdg->GetParticle(mctrack->
GetPdgCode())->Charge() << endl;
344 for (
int r = 0;
r < recotrackmcid.size();
r++) {
345 if (recotrackmcid.at(
r) == j) {
348 if (calcphimc < 0) calcphimc += 2*
TMath::Pi();
349 Double_t calcphireco = recotrackphi.at(
r);
350 if (calcphireco < 0) calcphireco += 2*
TMath::Pi();
351 Double_t deltaphi = calcphireco - calcphimc;
358 Double_t deltapt = recotrackpt.at(
r) - mcpt;
360 Int_t chargesign = 1;
361 if (pdg->GetParticle(mctrack->
GetPdgCode())->Charge() < 0) {
365 hphires->Fill(deltaphi*chargesign);
366 hptres->Fill(mcpt, deltapt/mcpt);
368 bool phiok = (deltaphi < 0.1);
373 bool ptok = ( (deltapt < mcpt * 0.1) || ((mcpt > 450)&&(recotrackpt.at(
r) > 450)) );
376 recotrackstatus.at(
r) = 1;
381 if ((!phiok) && ptok) {
384 if (phiok && (!ptok)) {
387 if ((!phiok) && (!ptok)) {
393 if (trackstatus == 0) {
399 mctrackstatus.push_back(trackstatus);
409 for (
int s = 0;
s < mctrackstatus.size();
s++) {
410 if (mctrackstatus.at(
s) > 0) {
415 int recotrackstatusgoodcount = 0;
416 int recotrackstatusbadcount = 0;
418 for (
int s = 0;
s < recotrackstatus.size();
s++) {
419 if (recotrackmcid.at(
s) <= globalmaxmcid) {
420 if (recotrackstatus.at(
s) > 0) {
421 recotrackstatusgoodcount++;
423 recotrackstatusbadcount++;
432 cout <<
"Total reconstructed track count: " << totalrecotracks << endl;
433 cout <<
"Reached MC ID: " << globalmaxmcid << endl;
434 cout <<
"Reconstructable MC track count: " << goodmctracks << endl;
435 cout <<
"MC Tracks reconstructed: " << goodrecotracks << endl;
436 cout <<
"MC Tracks NOT reconstructed: " << norecotracks << endl;
437 cout <<
"Reco Status Good: " << recotrackstatusgoodcount << endl;
438 cout <<
"Reco Status Bad: " << recotrackstatusbadcount << endl;
439 cout <<
"Delta phi failed: " << deltaphifailed << endl;
440 cout <<
"Delta pt failed: " << deltaptfailed << endl;
441 cout <<
"Delta both failed: " << deltabothfailed << endl;
442 cout <<
"Secondary Tracks: " << secondaries << endl;
444 cout <<
"MC Match Good 90: " << mcidgood90 << endl;
445 cout <<
"MC Match Good 80: " << mcidgood80 << endl;
446 cout <<
"MC Match Bad: " << mcidbad << endl;
450 TFile outputstorage(histofilename,
"UPDATE");
473 cout << endl << endl;
474 cout <<
"Macro finished succesfully." << endl;
475 cout <<
"Real time " << rtime <<
" s, CPU time " << ctime <<
" s" << endl;
TFile filereco("MvdStt_Test_reco.root")
Int_t GetNPoints(DetectorId detId) const
TVector3 GetMomentum() const
TFile filedigi("testdigi.root")
TVector3 GetStartVertex() const
Int_t GetMotherID() const
int runOnlineDisplayMCCheck(Int_t maximumTime=2050)
void SaveAndUpdateHisto(TH1 *currenthisto, TFile &storagefile)