FairRoot/PandaRoot
TB_checkHitReconstruction.C
Go to the documentation of this file.
1 int TB_checkHitReconstruction(Double_t pntHitDist=0.1) {
2 
3  TFile* simFile = TFile::Open("Gem_MvdStt_3Stations_DPM_n10000.root");
4  TTree* simTree = (TTree*)simFile->Get("pndsim");
5  TFile* digFile = TFile::Open("Gem_MvdStt_3Stations_DPM_n10000_digiSorted_a10000_c1000x10.root");
6  TTree* digTree = (TTree*)digFile->Get("pndsim");
7  TFile* cluFile = TFile::Open("Gem_MvdStt_3Stations_DPM_n10000_CluHi_a10000_c1000x10.root");
8  TTree* cluTree = (TTree*)cluFile->Get("pndsim");
9 
10  cout << "FILES FETCHED" << endl;
11 
12  FairEventHeader* evHeader = NULL;
13  TClonesArray* mcTrackArray = new TClonesArray("PndMCTrack");
14  TClonesArray* gemPointArray = new TClonesArray("PndGemMCPoint");
15  TClonesArray* gemHitArray = new TClonesArray("PndGemHit");
16 
17  simTree->SetBranchAddress("MCTrack",&mcTrackArray);
18  simTree->SetBranchAddress("GEMPoint",&gemPointArray);
19  digTree->SetBranchAddress("EventHeader.",&evHeader);
20  cluTree->SetBranchAddress("GEMHit",&gemHitArray);
21 
22  cout << "GOT ALL THE BRANCHES" << endl;
23 
24  TH2F* hist_YT = new TH2F("hist_YT","hist_YT",200,-10,10,200,-75,125);
25  TH2F* hist_XT = new TH2F("hist_XT","hist_XT",200,-10,10,200,-75,125);
26  TH2F* hist_XY = new TH2F("hist_XY","hist_XY",200,-10,10,200,-10,10);
27  TH1F* hist_T = new TH1F("hist_T", "hist_T", 500,-200,300);
28  TH1F* hist_Y = new TH1F("hist_Y", "hist_Y", 500,-25,25);
29  TH1F* hist_X = new TH1F("hist_X", "hist_X", 500,-25,25);
30  TH1F* hist_N = new TH1F("hist_N","Nof of matched hits",101,-0.5,100.5);
31 
32  const Int_t kNofStat = 3;
33 
34  const Int_t kNofEvents = 10000;//cluTree->GetEntries();
35  const Int_t kMaxNofHits = 10000;
36  // cout << "what's not static here?"<< endl;
37  Bool_t hitAssignedToPoint[kNofEvents][kMaxNofHits];
38  // cout << "really?" << endl
39  for ( Int_t iev = 0 ; iev < kNofEvents ; iev++ ) {
40  for ( Int_t ihit = 0 ; ihit < kMaxNofHits ; ihit++ ) {
41  hitAssignedToPoint[iev][ihit] = kFALSE;
42  }
43  }
44  cout << "DONE" << endl;
45 
46  Double_t statZPos[kNofStat] = {117.,153.,189.};
47 
48  TH1F* fhPointRadNof [kNofStat][2];
49  TH1F* fhPointRadReco [kNofStat][2];
50  TH1F* fhPointRadRecoEff[kNofStat][2];
51  TH2F* fhPointNof [kNofStat][2];
52  TH2F* fhPointReco [kNofStat][2];
53  TH2F* fhPointRecoEff [kNofStat][2];
54 
55  TH1F* fhHitRadNof [kNofStat][2];
56  TH1F* fhHitRadFake [kNofStat][2];
57  TH1F* fhHitRadFakeProb[kNofStat][2];
58  TH2F* fhHitNof [kNofStat][2];
59  TH2F* fhHitFake [kNofStat][2];
60  TH2F* fhHitFakeProb [kNofStat][2];
61 
62  for ( Int_t istat = 0 ; istat < kNofStat ; istat++ ) {
63  for ( Int_t isens = 0 ; isens < 2 ; isens++ ) {
64  fhPointRadNof [istat][isens] = new TH1F(Form("fhPointRadNof_s%d_s%d",istat,isens),
65  Form("Number of points, station %d, sensor %d",istat,isens),
66  100,0,100);
67  fhPointRadReco [istat][isens] = new TH1F(Form("fhPointRadReco_s%d_s%d",istat,isens),
68  Form("Number of reconstructed points (point-hit %.2f cm), station %d, sensor %d",pntHitDist,istat,isens),
69  100,0,100);
70  fhPointRadRecoEff [istat][isens] = new TH1F(Form("fhPointRadRecoEff_s%d_s%d",istat,isens),
71  Form("Hit finding efficiency (point-hit %.2f cm), station %d, sensor %d",pntHitDist,istat,isens),
72  100,0,100);
73  fhPointNof [istat][isens] = new TH2F(Form("fhPointNof_s%d_s%d",istat,isens),
74  Form("Number of points, station %d, sensor %d",istat,isens),
75  140,-70,70,
76  140,-70,70);
77  fhPointReco [istat][isens] = new TH2F(Form("fhPointReco_s%d_s%d",istat,isens),
78  Form("Number of reconstructed points (point-hit %.2f cm), station %d, sensor %d",pntHitDist,istat,isens),
79  140,-70,70,
80  140,-70,70);
81  fhPointRecoEff [istat][isens] = new TH2F(Form("fhPointRecoEff_s%d_s%d",istat,isens),
82  Form("Hit finding efficiency (point-hit %.2f cm), station %d, sensor %d",pntHitDist,istat,isens),
83  140,-70,70,
84  140,-70,70);
85  fhHitRadNof [istat][isens] = new TH1F(Form("fhHitRadNof_s%d_s%d",istat,isens),
86  Form("Number of hits, station %d, sensor %d",istat,isens),
87  100,0,100);
88  fhHitRadFake [istat][isens] = new TH1F(Form("fhHitRadFake_s%d_s%d",istat,isens),
89  Form("Number of fake hits (point-hit %.2f cm), station %d, sensor %d",pntHitDist,istat,isens),
90  100,0,100);
91  fhHitRadFakeProb [istat][isens] = new TH1F(Form("fhHitRadFakeProb_s%d_s%d",istat,isens),
92  Form("Fake hit probability (point-hit %.2f cm), station %d, sensor %d",pntHitDist,istat,isens),
93  100,0,100);
94  fhHitNof [istat][isens] = new TH2F(Form("fhHitNof_s%d_s%d",istat,isens),
95  Form("Number of hits, station %d, sensor %d",istat,isens),
96  140,-70,70,
97  140,-70,70);
98  fhHitFake [istat][isens] = new TH2F(Form("fhHitFake_s%d_s%d",istat,isens),
99  Form("Number of fake hits (point-hit %.2f cm), station %d, sensor %d",pntHitDist,istat,isens),
100  140,-70,70,
101  140,-70,70);
102  fhHitFakeProb [istat][isens] = new TH2F(Form("fhHitFakeProb_s%d_s%d",istat,isens),
103  Form("Fake hit probability (point-hit %.2f cm), station %d, sensor %d",pntHitDist,istat,isens),
104  140,-70,70,
105  140,-70,70);
106  }
107  }
108 
109  PndGemMCPoint* tempPoint = NULL;
110  PndGemHit* tempHit = NULL;
111 
112  Int_t startRev = 0;
113 
114  for ( Int_t imev = 0 ; imev < simTree->GetEntries() ; imev++ ) {
115  // for ( Int_t imev = 0 ; imev < 100 ; imev++ ) {
116  simTree->GetEntry(imev);
117  digTree->GetEntry(imev);
118  cout << "\r MC EVENT " << imev << flush;
119 
120  for ( Int_t ipoint = 0 ; ipoint < gemPointArray->GetEntries() ; ipoint++ ) {
121  tempPoint = (PndGemMCPoint*)gemPointArray->At(ipoint);
122  // cout << "POINT @ ( " << tempPoint->GetX() << " , " << tempPoint->GetY() << " , " << tempPoint->GetZ() << " )"
123  // << " at time " << evHeader->GetEventTime() << " + " << tempPoint->GetTime() << " ns" << endl;
124  Int_t pntStat = -1;
125  for ( Int_t istat = 0 ; istat < kNofStat ; istat++ ) {
126  if ( TMath::Abs(statZPos[istat]-tempPoint->GetZ()) < 5. ) {
127  pntStat = istat;
128  }
129  }
130  if ( pntStat == -1 ) cout << "Station at " << tempPoint->GetZ() << " not identified" << endl;
131 
132  Int_t pntSens = (tempPoint->GetZ()>=statZPos[pntStat]?1:0);
133  Double_t pntX = (tempPoint->GetX()+tempPoint->GetXOut())/2.;
134  Double_t pntY = (tempPoint->GetY()+tempPoint->GetYOut())/2.;
135  Double_t pntRad = TMath::Sqrt(pntX*pntX+pntY*pntY);
136 
137  Int_t irev = 0;
138 
139  // cout << "\r Event " << imev << " check from " << startRev << flush;
140  Int_t nofMatchedHits = 0;
141 
142  for ( irev = startRev ; irev < cluTree->GetEntries() ; irev++ ) {
143  // if ( verbose )
144  // if ( ipoint == 0 )
145  // cout << " EVENT " << irev << endl;
146 
147  Int_t nofUsableHits = 0;
148 
149  cluTree->GetEntry(irev);
150 
151  if ( gemHitArray->GetEntries() > 0 ) { // check if first hit has time much bigger than point time
152  tempHit = (PndGemHit*)gemHitArray->At(0);
153  if ( tempHit->GetTimeStamp() > evHeader->GetEventTime()+tempPoint->GetTime() + 200. )
154  break;
155  }
156 
157  if ( gemHitArray->GetEntries() > kMaxNofHits )
158  cout << "THERE ARE MORE HITS THAN EXPECTED, " << gemHitArray->GetEntries() << " IN EVENT " << irev << "." << endl;
159 
160  for ( Int_t ihit = 0 ; ihit < gemHitArray->GetEntries() ; ihit++ ) {
161  tempHit = (PndGemHit*)gemHitArray->At(ihit);
162 
163  if ( ipoint == 0 ) {
164  if ( tempHit->GetTimeStamp() > evHeader->GetEventTime() )
165  nofUsableHits ++;
166  }
167 
168  // cout << "HIT @ ( " << tempHit->GetX() << " , " << tempHit->GetY() << " , " << tempHit->GetZ() << " )"
169  // << " at time " << tempHit->GetTimeStamp() << endl;
170 
171  if ( TMath::Abs(tempHit->GetZ()-tempPoint->GetZ()) > 1. ) continue;// same station
172 
173  Bool_t goodT = kFALSE;
174  Bool_t goodX = kFALSE;
175  Bool_t goodY = kFALSE;
176 
177  // if ( TMath::Abs(tempHit->GetX()-tempPoint->GetX()) < pntHitDist ) goodX = kTRUE;
178  // if ( TMath::Abs(tempHit->GetY()-tempPoint->GetY()) < pntHitDist ) goodY = kTRUE;
179  if ( TMath::Abs(tempHit->GetX()-pntX) < pntHitDist ) goodX = kTRUE;
180  if ( TMath::Abs(tempHit->GetY()-pntY) < pntHitDist ) goodY = kTRUE;
181  if ( TMath::Abs(tempHit->GetTimeStamp()-(evHeader->GetEventTime()+tempPoint->GetTime())-5.5) < 8. ) goodT = kTRUE;
182 
183  if ( goodX ) {
184  hist_YT->Fill(tempHit->GetY()-pntY,
185  tempHit->GetTimeStamp()-(evHeader->GetEventTime()+tempPoint->GetTime()));
186  if ( goodY ) {
187  hist_T->Fill(tempHit->GetTimeStamp()-(evHeader->GetEventTime()+tempPoint->GetTime()));
188  if ( goodT ) {
189  // ((PndGemHit*)gemHitArray->At(ihit))->SetCharge(-100);
190  hitAssignedToPoint[irev][ihit] = kTRUE;
191  nofMatchedHits++;
192  }
193  }
194  if ( goodT )
195  hist_Y->Fill(tempHit->GetY()-pntY);//tempPoint->GetY());
196  }
197  if ( goodY ) {
198  hist_XT->Fill(tempHit->GetX()-pntX,//tempPoint->GetX(),
199  tempHit->GetTimeStamp()-(evHeader->GetEventTime()+tempPoint->GetTime()));
200  if ( goodT )
201  hist_X->Fill(tempHit->GetX()-pntX);//tempPoint->GetX());
202  }
203  if ( goodT ) {
204  hist_XY->Fill(tempHit->GetX()-pntX,//tempPoint->GetX(),
205  tempHit->GetY()-pntY);//tempPoint->GetY());
206  }
207 
208  }
209  if ( ipoint == 0 && nofUsableHits == 0 && gemHitArray->GetEntries() > 0 )
210  startRev++;
211  }
212  hist_N->Fill(nofMatchedHits);
213 
214  fhPointRadNof[pntStat][pntSens]->Fill(pntRad);
215  fhPointNof [pntStat][pntSens]->Fill(pntX,pntY);
216  if ( nofMatchedHits > 0 ) {
217  fhPointRadReco[pntStat][pntSens]->Fill(pntRad);
218  fhPointReco [pntStat][pntSens]->Fill(pntX,pntY);
219  }
220  }
221  }
222 
223  for ( irev = 0 ; irev < cluTree->GetEntries() ; irev++ ) {
224  // for ( irev = 0 ; irev < 100 ; irev++ ) {
225  cluTree->GetEntry(irev);
226  for ( Int_t ihit = 0 ; ihit < gemHitArray->GetEntries() ; ihit++ ) {
227  //cout << "event " << irev << " hit " << ihit << " - " << tempHit->GetCharge() << " ----> " << (hitAssignedToPoint[irev][ihit]?"MATCHED":"") << end;l
228  tempHit = (PndGemHit*)gemHitArray->At(ihit);
229 
230  fhHitRadNof [tempHit->GetStationNr()-1][tempHit->GetSensorNr()-1]->Fill(TMath::Sqrt(tempHit->GetX()*tempHit->GetX()+tempHit->GetY()*tempHit->GetY()));
231  fhHitNof [tempHit->GetStationNr()-1][tempHit->GetSensorNr()-1]->Fill(tempHit->GetX(),tempHit->GetY());
232  if ( !hitAssignedToPoint[irev][ihit] ) {
233  fhHitRadFake [tempHit->GetStationNr()-1][tempHit->GetSensorNr()-1]->Fill(TMath::Sqrt(tempHit->GetX()*tempHit->GetX()+tempHit->GetY()*tempHit->GetY()));
234  fhHitFake [tempHit->GetStationNr()-1][tempHit->GetSensorNr()-1]->Fill(tempHit->GetX(),tempHit->GetY());
235  }
236  }
237  }
238 
239  for ( Int_t istat = 0 ; istat < kNofStat ; istat++ ) {
240  for ( Int_t isens = 0 ; isens < 2 ; isens++ ) {
241  fhPointRadNof [istat][isens]->Sumw2();
242  fhPointRadReco [istat][isens]->Sumw2();
243  fhPointRadRecoEff [istat][isens]->Divide(fhPointRadReco [istat][isens],fhPointRadNof [istat][isens],1.,1.,"B");
244  fhPointRadRecoEff [istat][isens]->Scale(100.);
245  fhHitRadNof [istat][isens]->Sumw2();
246  fhHitRadFake [istat][isens]->Sumw2();
247  fhHitRadFakeProb [istat][isens]->Divide(fhHitRadFake [istat][isens],fhHitRadNof [istat][isens],1.,1.,"B");
248  fhHitRadFakeProb [istat][isens]->Scale(100.);
249 
250  fhPointNof [istat][isens]->Sumw2();
251  fhPointReco [istat][isens]->Sumw2();
252  fhPointRecoEff[istat][isens]->Divide(fhPointReco[istat][isens],fhPointNof[istat][isens],1.,1.,"B");
253  fhPointRecoEff[istat][isens]->Scale(100.);
254  fhHitNof [istat][isens]->Sumw2();
255  fhHitFake [istat][isens]->Sumw2();
256  fhHitFakeProb [istat][isens]->Divide(fhHitFake [istat][isens],fhHitNof [istat][isens],1.,1.,"B");
257  fhHitFakeProb [istat][isens]->Scale(100.);
258  }
259  }
260 
261  TCanvas* canv_XY = new TCanvas("canv_XY","canv_XY",10,10,1000,800);
262  hist_XY->DrawClone("colz");
263 
264  TCanvas* canv_XT = new TCanvas("canv_XT","canv_XT",10,10,1000,800);
265  hist_XT->DrawClone("colz");
266 
267  TCanvas* canv_YT = new TCanvas("canv_YT","canv_YT",10,10,1000,800);
268  hist_YT->DrawClone("colz");
269 
270  TCanvas* canv_X = new TCanvas("canv_X","canv_X",10,10,1000,800);
271  hist_X->DrawClone();
272 
273  TCanvas* canv_Y = new TCanvas("canv_Y","canv_Y",10,10,1000,800);
274  hist_Y->DrawClone();
275 
276  TCanvas* canv_T = new TCanvas("canv_T","canv_T",10,10,1000,800);
277  hist_T->DrawClone();
278 
279  TCanvas* canv_N = new TCanvas("canv_N","canv_N",10,10,1000,800);
280  hist_N->DrawClone();
281 
282  TCanvas* canv_E = new TCanvas("canv_E","canv_E",10,10,1000,800);
283  TString drawString = "";
284  for ( Int_t istat = 0 ; istat < kNofStat ; istat++ ) {
285  for ( Int_t isens = 0 ; isens < 2 ; isens++ ) {
286  fhPointRadRecoEff [istat][isens]->SetLineColor(2+istat+isens*4);
287  fhPointRadRecoEff [istat][isens]->SetLineWidth(2);
288  fhPointRadRecoEff [istat][isens]->DrawClone(drawString.Data());
289  drawString = "kSAME";
290  }
291  }
292 
293  TCanvas* canv_F = new TCanvas("canv_F","canv_F",10,10,1000,800);
294  TString drawString = "";
295  for ( Int_t istat = 0 ; istat < kNofStat ; istat++ ) {
296  for ( Int_t isens = 0 ; isens < 2 ; isens++ ) {
297  fhHitRadFakeProb [istat][isens]->SetLineColor(2+istat+isens*4);
298  fhHitRadFakeProb [istat][isens]->SetLineWidth(2);
299  fhHitRadFakeProb [istat][isens]->DrawClone(drawString.Data());
300  drawString = "kSAME";
301  }
302  }
303 
304  TFile* outFile = new TFile(Form("hit_reco_%.1f.root",pntHitDist),"RECREATE");
305 
306  hist_YT->Write();
307  hist_XT->Write();
308  hist_XY->Write();
309  hist_T ->Write();
310  hist_Y ->Write();
311  hist_X ->Write();
312  hist_N ->Write();
313 
314  for ( Int_t istat = 0 ; istat < kNofStat ; istat++ ) {
315  for ( Int_t isens = 0 ; isens < 2 ; isens++ ) {
316  fhPointRadNof [istat][isens]->Write();
317  fhPointRadReco [istat][isens]->Write();
318  fhPointRadRecoEff[istat][isens]->Write();
319  fhHitRadNof [istat][isens]->Write();
320  fhHitRadFake [istat][isens]->Write();
321  fhHitRadFakeProb [istat][isens]->Write();
322  fhPointNof [istat][isens]->Write();
323  fhPointReco [istat][isens]->Write();
324  fhPointRecoEff [istat][isens]->Write();
325  fhHitNof [istat][isens]->Write();
326  fhHitFake [istat][isens]->Write();
327  fhHitFakeProb [istat][isens]->Write();
328  }
329  }
330  outFile->Write();
331  outFile->Close();
332  return 0;
333 }
outFile
TString outFile
Definition: hit_dirc.C:17
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndGemHit::GetSensorNr
Int_t GetSensorNr() const
Definition: PndGemHit.h:83
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
simFile
TString simFile
Definition: bump_emc.C:11
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
TB_checkHitReconstruction
int TB_checkHitReconstruction(Double_t pntHitDist=0.1)
Definition: TB_checkHitReconstruction.C:1
PndGemHit::GetStationNr
Int_t GetStationNr() const
Definition: PndGemHit.h:81
PndGemMCPoint::GetYOut
Double_t GetYOut() const
Definition: PndGemMCPoint.h:84
PndGemMCPoint::GetXOut
Double_t GetXOut() const
Definition: PndGemMCPoint.h:83
Generated on Wed Apr 3 2019 10:03:02 for FairRoot/PandaRoot by   doxygen 1.8.5