FEDRA emulsion software from the OPERA Collaboration
Functions | Variables
testFakeDoubletBTReduction.C File Reference

Functions

void ReadRunFakeDoubletBGReductionPlots ()
 

Variables

Int_t CUTTYPE = 0
 
Int_t SELECTTYPE = 3
 

Function Documentation

◆ ReadRunFakeDoubletBGReductionPlots()

void ReadRunFakeDoubletBGReductionPlots ( )

AATTENTION CHECKING function "angle" of EdbSegP (which is asin really)...

29 {
30
31 // Quick runthrough to check data integrity and so on...
32 Bool_t QUICK = 0;
33 // Read PVR object either from a root file (0) or from lnk.def data proc (default)
34 Bool_t READTYPE = 1;
35
36 TObjArray* ali_array = new TObjArray();
37 // maximal number of root files:
38 // (to be adapted in case of special needs)
39 Int_t N_array=12;
40 if (CUTTYPE==2) N_array=3;
41
42
43 Int_t StartArray=0;
44
45 EdbPVRec* ali;
46
47 if (READTYPE==0) {
48 for (Int_t n_array = StartArray; n_array < N_array; ++n_array) {
49 cout << "n_array = " << n_array << endl;
50 // This time, read from the pre-filtered file:
51 TFile* f = new TFile(Form("ScanVolume_Ali_%d_CUTTYPE_%d_SELECTTYPE_%d.root",n_array,CUTTYPE,999));
52 ali = new EdbPVRec();
53 ali = (EdbPVRec*)f->Get("EdbPVRec");
54 ali->Print();
55 ali_array->Add(ali);
56 }
57 } else {
58 ali = new EdbPVRec();
59 EdbDataProc* proc = new EdbDataProc("lnk.def");
60 proc->InitVolume(ali,0);
61 cout << "read_pvr.C :: I have read a EdbPVRec* object from basetrack cp data. " << endl;
62 cout << "read_pvr.C :: ali = " << ali << " and it has ali->Npatterns()= " << ali->Npatterns() << endl;
63 ali_array->Add(ali);
64 }
65
66//cout << "return HERE" << endl;
67//return;
68
69 // All Segment Pairs
70 TH2D* h_dR_dT = new TH2D("h_dR_dT","",200,0,0.2,200,0,40);
71 TH2D* h2_dR_dT = new TH2D("h2_dR_dT","",200,0,0.2,200,0,40);
72 // after the cut into the line area (to see if these BTs are identical,
73 // so should also be their chi2 and W values)s
74 TH2D* h_dchi2_dW = new TH2D("h_dchi2_dW","",13,-6.5,6.5,100,-1,1);
75 TH2D* h2_dchi2_dW = new TH2D("h2_dchi2_dW","",13,-6.5,6.5,100,-1,1);
76
77 TH2D* h_chi2_W = new TH2D("h_chi2_W","",26,5.5,31.5,100,0,2);
78 TH2D* h2_chi2_W = new TH2D("h2_chi2_W","",26,5.5,31.5,100,0,2);
79
80 TH2D* h_TX_TY = new TH2D("h_TX_TY","",100,-0.5,0.5,100,-0.5,0.5);
81 TH2D* h2_TX_TY = new TH2D("h2_TX_TY","",100,-0.5,0.5,100,-0.5,0.5);
82
83 TH2D* h_DeltaTanTheta_Angle = new TH2D("h_DeltaTanTheta_Angle","",200,0,0.2,200,0,0.2);
84
85 TH2D* h_X_Y = new TH2D("h_X_Y","",1250,0,125000,1000,0,100000);
86 TH2D* h2_X_Y = new TH2D("h2_X_Y","",1250,0,125000,1000,0,100000);
87
88 TH2D* h_X_Y = new TH2D("h_X_Y","",2*1250,0,125000,2*1000,0,100000);
89 TH2D* h2_X_Y = new TH2D("h2_X_Y","",2*1250,0,125000,2*1000,0,100000);
90
91 TH1D* h2_X_Y_ProjX = new TH1D("h2_X_Y_ProjX","",2*1250,0,125000);
92
93
94 TCanvas* c_dR_dT = new TCanvas();
95 TCanvas* c_dchi2_dW = new TCanvas();
96 TCanvas* c2_dR_dT = new TCanvas();
97 TCanvas* c2_dchi2_dW = new TCanvas();
98 TCanvas* c_chi2_W = new TCanvas();
99 TCanvas* c2_chi2_W = new TCanvas();
100 TCanvas* c_TX_TY = new TCanvas("c_TX_TY","c_TX_TY",2);
101 TCanvas* c2_TX_TY = new TCanvas("c2_TX_TY","c2_TX_TY",2);
102 TCanvas* c_DeltaTanTheta_Angle = new TCanvas("c_DeltaTanTheta_Angle","c_DeltaTanTheta_Angle");
103
104 Float_t deltaR=0;
105 Float_t deltaTheta=0;
106 Float_t deltachi2=0;
107 Float_t deltaW=0;
108 Float_t deltaTanTheta=0;
109 Float_t angleTheta=0;
110
111 // Bool Variables for Region
112 Bool_t IsRegion0;
113 Bool_t IsRegion1;
114 Bool_t IsRegion2;
115
116 // Create Variables For ExtractSubpattern boundaries
117 Float_t mini[5];
118 Float_t maxi[5];
119 mini[2]=-1;
120 mini[3]=-1;
121 mini[4]=0.0;
122 maxi[2]=1;
123 maxi[3]=1;
124 maxi[4]=32.0;
125 Float_t halfpatternsize=40;
126
127 Double_t deltaThetaTEST;
128 Double_t dR_allowed;
129 Double_t dR_allowed_m_deltaR;
130
131 int n_array = 0;
132
133 TCanvas* c_X_Y = new TCanvas(Form("c_X_Y_%d",n_array),Form("c_X_Y_%d",n_array));
134 TCanvas* c2_X_Y = new TCanvas(Form("c2_X_Y_%d",n_array),Form("c2_X_Y_%d",n_array));
135
136
137 for (n_array = 0; n_array < N_array; ++n_array) {
138// for (n_array = 0; n_array < 1; ++n_array) {
139
140 cout << "Doing Array " << n_array << " ----------------------------------------" << endl;
141 EdbPVRec* ali = (EdbPVRec*)ali_array->At(n_array);
142 aliSource = ali;
143
145 aliNEW = (EdbPVRec*)ali->Clone();
146 aliNEW->Print();
147 Int_t npataliNEW = aliNEW->Npatterns();
148 for (int ctr_i = 0; ctr_i< npataliNEW; ++ctr_i ) aliNEW->GetPattern(ctr_i)->Reset();
150
151 /*
152 h_X_Y->Reset();
153 h2_X_Y->Reset();
154 TCanvas* c_X_Y = new TCanvas(Form("c_X_Y_%d",n_array),Form("c_X_Y_%d",n_array));
155 TCanvas* c2_X_Y = new TCanvas(Form("c2_X_Y_%d",n_array),Form("c2_X_Y_%d",n_array));
156 */
157
158 for (int i = 0; i <aliSource->Npatterns(); i++ ) {
159
160 cout << "Doing Pattern " << i << endl;
161
162 EdbPattern* pat = aliSource->GetPattern(i);
163 Int_t nPat=pat->N();
164
165 EdbPattern* patClone = aliNEW->GetPattern(i);
166 Int_t nPatClone=patClone->N();
167
168 for (int j = 0; j < nPat; j++ ) {
169
170 if (j%25000==0) cout << j << " out of " << nPat << endl;
171
172 seg = pat->GetSegment(j);
173// if (seg->MCEvt()>0) continue;
174
175 mini[0]=seg->X()-halfpatternsize;
176 mini[1]=seg->Y()-halfpatternsize;
177 maxi[0]=seg->X()+halfpatternsize;
178 maxi[1]=seg->Y()+halfpatternsize;
179 EdbPattern* singlePattern=(EdbPattern*)aliSource->GetPattern(i)->ExtractSubPattern(mini,maxi,0);
180 Int_t nPatsinglePattern=singlePattern->N();
181 /*
182 cout <<"------------------------"<< endl;
183 cout << j << " out of " << nPat << endl;
184 cout << nPatsinglePattern << " nPatsinglePattern entries for this subpattern" << endl;
185 singlePattern->Print();
186 cout <<"------------------------"<< endl;
187 */
188 for (int k = 0; k < nPatsinglePattern; k++ ) {
189
190 EdbSegP* seg1 = singlePattern->GetSegment(k);
191
192// if (seg1->MCEvt()>0) continue;
193 // if (seg==seg1) continue; // doesnt work, since the new subpattern has new adresses
194 // for the segments, so rather compare segments with the function EdbSegP::IsEqual
195 if (seg->IsEqual(seg1)) continue;
196
197 // Skip already large distances
198 if (TMath::Abs(seg->X()-seg1->X())>40.1) continue;
199 if (TMath::Abs(seg->Y()-seg1->Y())>40.1) continue;
200 if (TMath::Abs(seg->TX()-seg1->TX())>0.21) continue;
201 if (TMath::Abs(seg->TY()-seg1->TY())>0.21) continue;
202
203 // To avoid double filling (seg-seg1 and seg1-seg)
204 // we restrict to seg->ID() < seg1->ID()
205 if (seg->ID() > seg1->ID()) continue;
206
207
208 // Take only interesting combinations:
209 deltaR=TMath::Sqrt( TMath::Power(seg->X()-seg1->X(),2)+TMath::Power(seg->Y()-seg1->Y(),2) );
210 if (deltaR>40) continue;
211
212 deltaTanTheta=TMath::Sqrt( TMath::Power(seg->TX()-seg1->TX(),2)+TMath::Power(seg->TY()-seg1->TY(),2) );
213 if (deltaTanTheta>0.2) continue;
214
216 deltaThetaTEST = TMath::ASin(seg->Angle(*seg,*seg1));
217
218
219 angleTheta = EdbMath::Angle3( seg->TX(), seg->TY(),seg1->TX(), seg1->TY() );
220 deltaTheta = deltaTanTheta;
221 // we explicitly use here not the (correct) angle,
222 // but the Delta Tan Version, as also in the shower reconstruction, to treat it in a consistent way!
223 deltachi2 = seg->Chi2()-seg1->Chi2();
224 deltaW = seg->W()-seg1->W();
225
226
227 // Fill Histograms (only for first segment,
228 // where value of just one segment is needed)
229 h_dR_dT->Fill(deltaTheta,deltaR);
230 h_dchi2_dW->Fill(deltaW,deltachi2);
231 h_chi2_W->Fill(seg->W(),seg->Chi2());
232 h_TX_TY->Fill(seg->TX(),seg->TY());
233 h_DeltaTanTheta_Angle ->Fill(angleTheta,deltaTanTheta);
234 h_X_Y->Fill(seg->X(),seg->Y());
235
236
237
238 // Region check of the (Fake-)BTPair
239 IsRegion0=kFALSE;
240 IsRegion1=kFALSE;
241 IsRegion2=kFALSE;
242
243 // Linear Correlation Line by experimental values:
244 // dR_allowed = 11(micron)/0.1(rad) * dTheta // b??????
245 // dR_allowed = 5(micron)/0.05(rad) * dTheta // b138756
246// dR_allowed = 5/0.05*deltaTheta;
247// dR_allowed = 10.75/0.1*deltaTheta;
248 dR_allowed = 14.0/0.13*deltaTheta;
249// 14(micron)/0.13(rad)
250 dR_allowed_m_deltaR=TMath::Abs(deltaR-dR_allowed);
251
252 // Set flags for region
253 if (dR_allowed_m_deltaR<0.5) IsRegion0=kTRUE;
254 if (deltaR<3.0 && deltaTheta<0.01) IsRegion1=kTRUE;
255 if (abs(deltaR-6.0)<2 && deltaTheta<0.01) IsRegion2=kTRUE;
256
257
258
259 //==================================================
260 if (SELECTTYPE==999) ; // dont impose extra cut...
261
262 if (SELECTTYPE==0) if (!IsRegion0) continue;
263
264 if (SELECTTYPE==1) if (!IsRegion1) continue;
265
266 if (SELECTTYPE==2) if (!IsRegion2) continue;
267
268 if (SELECTTYPE==3) if ( !(IsRegion0 || IsRegion1 || IsRegion2) ) continue;
269 //==================================================
270
271
272 // Fill Histograms (only for first segment,
273 // where value of just one segment is needed)
274 h2_dR_dT->Fill(deltaTheta,deltaR);
275 h2_dchi2_dW->Fill(deltaW,deltachi2);
276 h2_chi2_W->Fill(seg->W(),seg->Chi2());
277 h2_TX_TY->Fill(seg->TX(),seg->TY());
278 h2_X_Y->Fill(seg->X(),seg->Y());
279
280
281 //cout << " For j = " << j << " we have a segment pair, which is close! k=" << k << endl;
282 //seg->PrintNice();
283 //seg1->PrintNice();
284
285 //==================================================
286 patClone->AddSegmentNoDuplicate(*seg); // Add this first.
287 patClone->AddSegmentNoDuplicate(*seg1);
288 //==================================================
289
290
291 } // for (int k = j+1; k <pat->N(); k++ )
292
293
294 delete singlePattern;
295
296 }// for (int j)
297
298
299 if (QUICK) i+=aliSource->Npatterns()-1;
300 } // for (int i)
301
302
303// // // // // // // // // // // // // // // // // // // // // // // // // // // //
304 TFile* f = new TFile(Form("ScanVolume_Ali_%d_CUTTYPE_%d_SELECTTYPE_%d.root",n_array,CUTTYPE,SELECTTYPE),"RECREATE");
305 aliNEW -> Write();
306 aliNEW->Print();
307// // // // // // // // // // // // // // // // // // // // // // // // // // // //
308
309
310 } // for (int n_array = 0; n_array < N_array; ++n_array)
311
312
313
314
315
316 // Draw histograms
317 if (1==1) {
318 c_dR_dT->cd();
319 h_dR_dT->Draw("colz");
320 h_dR_dT->GetXaxis()->SetTitle("#Delta tan #theta");
321 h_dR_dT->GetYaxis()->SetTitle("#DeltaR (#mum)");
322 c_dR_dT->SetLogz();
323 gPad->Update();
324 c_dR_dT->Print("Plot_ViewOverlap_CPData_dR_dT.pdf");
325
326 c_dchi2_dW->cd();
327 h_dchi2_dW->Draw("colz");
328 gPad->Update();
329 h_dchi2_dW->GetYaxis()->SetTitle("#Delta#chi^{2}");
330 h_dchi2_dW->GetXaxis()->SetTitle("#DeltaW");
331 gPad->Update();
332 c_dchi2_dW->Print("Plot_ViewOverlap_CPData_dchi2_dW.pdf");
333
334 c2_dR_dT->cd();
335 h2_dR_dT->Draw("colz");
336 gPad->Update();
337 h2_dR_dT->GetXaxis()->SetTitle("#Delta tan #theta");
338 h2_dR_dT->GetYaxis()->SetTitle("#DeltaR (#mum)");
339 c2_dR_dT->SetLogz();
340 gPad->Update();
341 c2_dR_dT->Print(Form("Plot_ViewOverlap_CPData_dR_dT_SELECTTYPE_%d.pdf",SELECTTYPE));
342
343 c2_dchi2_dW->cd();
344 h2_dchi2_dW->Draw("colz");
345 h2_dchi2_dW->GetYaxis()->SetTitle("#Delta#chi^{2}");
346 h2_dchi2_dW->GetXaxis()->SetTitle("#DeltaW");
347 gPad->Update();
348 c2_dchi2_dW->Print(Form("Plot_ViewOverlap_CPData_dchi2_dW_SELECTTYPE_%d.pdf",SELECTTYPE));
349
350 c_chi2_W->cd();
351 h_chi2_W->Draw("colz");
352 h_chi2_W->GetYaxis()->SetTitle("#chi^{2}");
353 h_chi2_W->GetXaxis()->SetTitle("W");
354 gPad->Update();
355 c_chi2_W->Print("Plot_ViewOverlap_CPData_chi2_W.pdf");
356
357 c2_chi2_W->cd();
358 h2_chi2_W->Draw("colz");
359 h2_chi2_W->GetYaxis()->SetTitle("#chi^{2}");
360 h2_chi2_W->GetXaxis()->SetTitle("W");
361 gPad->Update();
362 c2_chi2_W->Print(Form("Plot_ViewOverlap_CPData_chi2_W_SELECTTYPE_%d.pdf",SELECTTYPE));
363
364 c_TX_TY->cd();
365 h_TX_TY->Draw("colz");
366 h_TX_TY->GetYaxis()->SetTitle("tan #theta_{Y}");
367 h_TX_TY->GetXaxis()->SetTitle("tan #theta_{X}");
368 c_TX_TY->SetLogz();
369 gPad->Update();
370 c_TX_TY->Print("Plot_ViewOverlap_CPData_TX_TY.pdf");
371
372 c2_TX_TY->cd();
373 h2_TX_TY->Draw("colz");
374 h2_TX_TY->GetYaxis()->SetTitle("tan #theta_{Y}");
375 h2_TX_TY->GetXaxis()->SetTitle("tan #theta_{X}");
376 c2_TX_TY->SetLogz();
377 gPad->Update();
378 c2_TX_TY->Print(Form("Plot_ViewOverlap_CPData_TX_TY_SELECTTYPE_%d.pdf",SELECTTYPE));
379
380
381 c_DeltaTanTheta_Angle->cd();
382 h_DeltaTanTheta_Angle->Draw("colz");
383 h_DeltaTanTheta_Angle->GetYaxis()->SetTitle("#Delta tan #theta");
384 h_DeltaTanTheta_Angle->GetXaxis()->SetTitle("angle #Theta (rad)");
385 c_DeltaTanTheta_Angle->SetLogz();
386 gPad->Update();
387 c_DeltaTanTheta_Angle->Print("Plot_ViewOverlap_CPData_DeltaTanTheta_Angle.pdf");
388
389
390 c_X_Y->cd();
391 h_X_Y->Draw("colz");
392 h_X_Y->GetXaxis()->SetTitle("X (#mum)");
393 h_X_Y->GetYaxis()->SetTitle("Y (#mum)");
394 c_X_Y->SetLogz();
395 gPad->Update();
396 TPaletteAxis *palette = (TPaletteAxis*)h_X_Y->GetListOfFunctions()->FindObject("palette");
397 palette->SetY1NDC(0.2);
398 palette->SetX2NDC(0.93);
399 gPad->Update();
400 c_X_Y->Print("Plot_ViewOverlap_CPData_X_Y.pdf");
401
402 c2_X_Y->cd();
403 h2_X_Y->Draw("colz");
404 h2_X_Y->GetXaxis()->SetTitle("X (#mum)");
405 h2_X_Y->GetYaxis()->SetTitle("Y (#mum)");
406 c2_X_Y->SetLogz();
407 gPad->Update();
408 TPaletteAxis *palette = (TPaletteAxis*)h2_X_Y->GetListOfFunctions()->FindObject("palette");
409 palette->SetY1NDC(0.2);
410 palette->SetX2NDC(0.93);
411 gPad->Update();
412 c2_X_Y->Print(Form("Plot_ViewOverlap_CPData_X_Y_SELECTTYPE_%d.pdf",SELECTTYPE));
413
414 } // Redraw Histos
415
416
417
418 TFile* f = new TFile(Form("HistosFakeDoubleBT_CUTTYPE_%d_SELECTTYPE_%d.root",CUTTYPE,SELECTTYPE),"RECREATE");
419 h_dR_dT -> Write();
420 h2_dR_dT-> Write();
421 h_dchi2_dW -> Write();
422 h2_dchi2_dW -> Write();
423 h_chi2_W -> Write();
424 h2_chi2_W -> Write();
425 h_TX_TY -> Write();
426 h2_TX_TY-> Write();
427 h_DeltaTanTheta_Angle -> Write();
428 h_X_Y -> Write();
429 h2_X_Y -> Write();
430 f->Close();
431
432 return;
433}
f
FILE * f
Definition: RecDispMC.C:150
EdbDataProc
emulsion data processing
Definition: EdbDataSet.h:181
EdbDataProc::InitVolume
int InitVolume(int datatype=0, const char *rcut="1")
Definition: EdbDataSet.cxx:2066
EdbMath::Angle3
static double Angle3(float tx1, float ty1, float tx2, float ty2)
Definition: EdbMath.cxx:42
EdbPVRec
Definition: EdbPVRec.h:148
EdbPattern
Definition: EdbPattern.h:273
EdbPattern::ExtractSubPattern
EdbPattern * ExtractSubPattern(float min[5], float max[5], int MCevt=-1)
Definition: EdbPattern.cxx:1470
EdbPatternsVolume::Npatterns
Int_t Npatterns() const
Definition: EdbPattern.h:366
EdbPatternsVolume::Print
void Print() const
Definition: EdbPattern.cxx:1693
EdbSegP
Definition: EdbSegP.h:21
EdbSegP::TX
Float_t TX() const
tangens = deltaX/deltaZ
Definition: EdbSegP.h:175
EdbSegP::ID
Int_t ID() const
Definition: EdbSegP.h:147
EdbSegP::X
Float_t X() const
Definition: EdbSegP.h:173
EdbSegP::Chi2
Float_t Chi2() const
Definition: EdbSegP.h:157
EdbSegP::Y
Float_t Y() const
Definition: EdbSegP.h:174
EdbSegP::W
Float_t W() const
Definition: EdbSegP.h:151
EdbSegP::TY
Float_t TY() const
tangens = deltaY/deltaZ
Definition: EdbSegP.h:176
EdbSegmentsBox::AddSegmentNoDuplicate
EdbSegP * AddSegmentNoDuplicate(EdbSegP &s)
Definition: EdbPattern.cxx:92
EdbSegmentsBox::N
Int_t N() const
Definition: EdbPattern.h:86
EdbSegmentsBox::GetSegment
EdbSegP * GetSegment(int i) const
Definition: EdbPattern.h:66
proc
EdbDataProc * proc
Definition: read_pvr.C:3
ali
EdbPVRec * ali
Definition: test_oracle.C:9
Write
ali Write()
TCanvas
new TCanvas()
SELECTTYPE
Int_t SELECTTYPE
Definition: testFakeDoubletBTReduction.C:25
CUTTYPE
Int_t CUTTYPE
Definition: testFakeDoubletBTReduction.C:24

Variable Documentation

◆ CUTTYPE

Int_t CUTTYPE = 0

◆ SELECTTYPE

Int_t SELECTTYPE = 3