CbmRoot
CbmRichGeoOpt.cxx
Go to the documentation of this file.
1 
8 #include "CbmRichGeoOpt.h"
9 #include "TCanvas.h"
10 #include "TClonesArray.h"
11 #include "TH1.h"
12 #include "TH1D.h"
13 #include "TH2D.h"
14 #include "TH3D.h"
15 
16 #include "CbmDrawHist.h"
17 #include "CbmMCTrack.h"
18 #include "CbmRichGeoManager.h"
19 #include "CbmRichHit.h"
20 #include "CbmRichRing.h"
21 #include "CbmTrackMatchNew.h"
22 #include "FairTrackParam.h"
23 
24 #include "CbmRichRingLight.h"
25 #include <boost/assign/list_of.hpp>
26 #include <iostream>
27 #include <string>
28 
29 using namespace std;
30 using boost::assign::list_of;
31 
32 CbmRichGeoOpt::CbmRichGeoOpt()
33  : FairTask("CbmRichGeoOpt")
34  , fRichPoints(NULL)
35  , fMcTracks(NULL)
36  , fRefPoints(NULL)
37  , fRichHits(NULL)
38  , fRichRings(NULL)
39  , fRichRingMatches(NULL)
40  , fEventNum(0)
41  , PMTPointsFilled(0)
42  , fMinNofHits(7)
43  , nPhotonsNotOnPlane(0)
44  , nPhotonsNotOnSphere(0)
45  , nTotalPhorons(0)
46  , PMTPlanePoints()
47  , PMTPlaneCenter()
48  , ReadPMTPlaneCenter()
49  , ReadPMTPlaneCenterOrig()
50  , MirrorCenter()
51  ,
52  // ReadMirrorCenter(),
53  RotX(0.)
54  , RotY(0.)
55  , PMT_r1()
56  , PMT_r2()
57  , PMT_norm()
58  , PMTPlaneCenterX(0.)
59  , PMTPlaneCenterY(0.)
60  , PMTPlaneThirdX(0.)
61  , MirrPosition()
62  , SensPointsFilled(0)
63  , SensPlanePoints()
64  , Sens_r1()
65  , Sens_r2()
66  , Sens_norm()
67  , H_Diff_LineRefPMT_MomAtPMT(NULL)
68  , H_Theta_TwoVectors(NULL)
69  , H_DistancePMTtoMirrCenter(NULL)
70  , H_DistancePMTtoMirr(NULL)
71  , H_MomPrim(NULL)
72  , H_PtPrim(NULL)
73  , H_MomPt(NULL)
74  , H_Mom_Theta_MC(NULL)
75  , H_Pt_Theta_MC(NULL)
76  , H_Mom_Theta_Rec(NULL)
77  , H_Pt_Theta_Rec(NULL)
78  , H_Mom_Theta_Acc(NULL)
79  , H_Pt_Theta_Acc(NULL)
80  , H_MomRing(NULL)
81  , H_acc_mom_el(NULL)
82  , H_acc_pty_el(NULL)
83  , H_Mom_XY_Theta25(NULL)
84  , H_MomPrim_RegularTheta(NULL)
85  , H_acc_mom_el_RegularTheta(NULL)
86  , H_Hits_XY(NULL)
87  , H_Hits_XY_LeftHalf(NULL)
88  , H_Hits_XY_RightHalf(NULL)
89  , H_Hits_XY_Left2Thirds(NULL)
90  , H_Hits_XY_RightThird(NULL)
91  , H_PointsIn_XY(NULL)
92  , H_PointsIn_XY_LeftHalf(NULL)
93  , H_PointsIn_XY_RightHalf(NULL)
94  , H_PointsIn_XY_Left2Thirds(NULL)
95  , H_PointsIn_XY_RightThird(NULL)
96  , H_PointsOut_XY(NULL)
97  , H_NofPhotonsPerEv(NULL)
98  , H_NofPhotonsPerHit(NULL)
99  , H_NofPhotonsSmallerThan30(NULL)
100  , H_DiffXhit(NULL)
101  , H_DiffYhit(NULL)
102  , H_dFocalPoint_Delta(NULL)
103  , H_dFocalPoint_Rho(NULL)
104  , H_Alpha(NULL)
105  , H_Alpha_UpLeft(NULL)
106  , H_Alpha_UpLeft_II(NULL)
107  , H_Alpha_UpLeft_RegularTheta(NULL)
108  , H_Alpha_UpLeft_LeftHalf(NULL)
109  , H_Alpha_UpLeft_RightHalf(NULL)
110  , H_Alpha_UpLeft_RightThird(NULL)
111  , H_Alpha_UpLeft_Left2Thirds(NULL)
112  , H_Alpha_XYposAtDet(NULL)
113  , H_Alpha_XYposAtDet_RegularTheta(NULL)
114  , H_Alpha_XYposAtDet_LeftHalf(NULL)
115  , H_Alpha_XYposAtDet_RightHalf(NULL)
116  , H_Alpha_XYposAtDet_RightThird(NULL)
117  , H_Alpha_XYposAtDet_Left2Thirds(NULL)
118  , H_NofHitsAll(NULL)
119  , H_NofRings(NULL)
120  , H_NofRings_NofHits(NULL)
121  , H_RingCenterX(NULL)
122  , H_RingCenterY(NULL)
123  , H_RingCenter(NULL)
124  , H_Radius(NULL)
125  , H_aAxis(NULL)
126  , H_bAxis(NULL)
127  , H_boa(NULL)
128  , H_boa_RegularTheta(NULL)
129  , H_boa_LeftHalf(NULL)
130  , H_boa_RightHalf(NULL)
131  , H_boa_RightThird(NULL)
132  , H_boa_Left2Thirds(NULL)
133  , H_dR(NULL)
134  , H_dR_aa(NULL)
135  , H_dR_RegularTheta(NULL)
136  , H_dR_LeftHalf(NULL)
137  , H_dR_RightHalf(NULL)
138  , H_dR_RightThird(NULL)
139  , H_dR_Left2Thirds(NULL)
140  , H_RingCenter_Aaxis(NULL)
141  , H_RingCenter_Baxis(NULL)
142  , H_RingCenter_boa(NULL)
143  , H_RingCenter_boa_RegularTheta(NULL)
144  , H_RingCenter_boa_LeftHalf(NULL)
145  , H_RingCenter_boa_RightHalf(NULL)
146  , H_RingCenter_boa_RightThird(NULL)
147  , H_RingCenter_boa_Left2Thirds(NULL)
148  , H_RingCenter_dR(NULL)
149  , H_RingCenter_dR_RegularTheta(NULL)
150  , H_RingCenter_dR_LeftHalf(NULL)
151  , H_RingCenter_dR_RightHalf(NULL)
152  , H_RingCenter_dR_RightThird(NULL)
153  , H_RingCenter_dR_Left2Thirds(NULL)
154 
155 {
156  /*
157  fEventNum = 0;
158  PMTPointsFilled = 0;
159  fMinNofHits = 7;
160  nPhotonsNotOnPlane = 0;
161  nPhotonsNotOnSphere = 0;
162  nTotalPhorons = 0;
163  */
164 }
165 
166 CbmRichGeoOpt::~CbmRichGeoOpt() {}
167 
168 InitStatus CbmRichGeoOpt::Init() {
169  cout << "CbmRichGeoOpt::Init" << endl;
170  FairRootManager* ioman = FairRootManager::Instance();
171  if (NULL == ioman) {
172  Fatal("CbmRichGeoOpt::Init", "RootManager not instantised!");
173  }
174 
175  fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
176  if (NULL == fRichPoints) {
177  Fatal("CbmRichGeoOpt::Init", "No RichPoint array!");
178  }
179 
180  fRefPoints = (TClonesArray*) ioman->GetObject("RefPlanePoint");
181  if (NULL == fRefPoints) {
182  Fatal("CbmRichGeoOpt::Init", "No fRefPoints array!");
183  }
184 
185  fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
186  if (NULL == fMcTracks) { Fatal("CbmRichGeoOpt::Init", "No MCTrack array!"); }
187 
188  fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
189  if (NULL == fRichHits) { Fatal("CbmRichGeoTest::Init", "No RichHit array!"); }
190 
192  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
193  if (NULL == fRichRings) {
194  Fatal("CbmRichGeoTest::Init", "No RichRing array!");
195  }
196 
197  fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
198  if (NULL == fRichRingMatches) {
199  Fatal("CbmRichGeoTest::Init", "No RichRingMatch array!");
200  }
201 
203  //cout<<" initializing points values -1000"<<endl;
204  PMTPlanePoints.resize(3);
205  SensPlanePoints.resize(3);
206  for (unsigned int p = 0; p < PMTPlanePoints.size(); p++) {
207  PMTPlanePoints[p].SetX(-1000.);
208  PMTPlanePoints[p].SetY(-1000.);
209  PMTPlanePoints[p].SetZ(-1000.);
210  SensPlanePoints[p].SetX(-1000.);
211  SensPlanePoints[p].SetY(-1000.);
212  SensPlanePoints[p].SetZ(-1000.);
213  }
214  // CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
215  PMTPlaneCenterX = 0.; //-1.*gp->fPmtWidthX/2.;
216  PMTPlaneCenterY = 0.; //-1.*gp->fPmtWidthY/2.;
217  PMTPlaneThirdX = 0.; //-1.*gp->fPmtWidthX/3.;
218  //cout<<"fGP.fPmtWidthX = "<< gp->fPmtWidthX<<", PMTPlaneCenterX = "<< PMTPlaneCenterX<<", PMTPlaneThirdX = "<< PMTPlaneThirdX<<endl;
219 
220  InitHistograms();
221  return kSUCCESS;
222 }
223 
224 void CbmRichGeoOpt::Exec(Option_t* /*option*/) {
225  fEventNum++;
226  // cout << "#################### CbmRichGeoOpt, event No. " << fEventNum << endl;
227 
228  if (PMTPointsFilled == 0) {
229  for (unsigned int p = 1; p < PMTPlanePoints.size(); p++) {
230  if (PMTPlanePoints[p].X() == PMTPlanePoints[p - 1].X()) {
231  cout << "PMTPlanePoints[" << p << "].X == PMTPlanePoints[" << p - 1
232  << "].X ==" << PMTPlanePoints[p - 1].X() << endl;
233  cout << "Calling FillPointsAtPMT()" << endl;
234  FillPointsAtPMT();
235  PMTPointsFilled = 0;
236  } else {
237  PMTPointsFilled = 1;
238  }
239  }
240  }
241  if (SensPointsFilled == 0) {
242  for (unsigned int p = 1; p < SensPlanePoints.size(); p++) {
243  if (SensPlanePoints[p].X() == SensPlanePoints[p - 1].X()) {
244  FillPointsAtPMTSensPlane();
245  SensPointsFilled = 0;
246  } else {
247  SensPointsFilled = 1;
248  }
249  }
250  }
251 
252  //cout << "#################### CbmRichGeoOpt, event No. " << fEventNum << endl;
253  //Fill the coordinates of the three points on the PMT plane
254  CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
255  PMTPlaneCenter.SetX(gp->fPmt.fX);
256  PMTPlaneCenter.SetY(gp->fPmt.fY);
257  PMTPlaneCenter.SetZ(gp->fPmt.fZ);
259  // if(PMTPointsFilled==1 && SensPointsFilled==1 ){cout<<" Both are filled."<<endl;}
260  // else if(PMTPointsFilled==1 && SensPointsFilled==0){cout<<" PMTPointsFilled==1 && SensPointsFilled==0."<<endl;}
261  // else if(PMTPointsFilled==0 && SensPointsFilled==1){cout<<" PMTPointsFilled==0 && SensPointsFilled==1."<<endl;}
262  // else{cout<<" PMTPointsFilled==0 && SensPointsFilled==0."<<endl;}
263 
264  if (PMTPointsFilled == 1) {
265  if (fEventNum < 10) {
266  for (unsigned int p = 0; p < PMTPlanePoints.size(); p++) {
267  cout << "Point " << p << ": (" << PMTPlanePoints[p].X() << " , "
268  << PMTPlanePoints[p].Y() << " , " << PMTPlanePoints[p].Z() << ")"
269  << endl;
270  }
271 
272  PMT_r1 = PMTPlanePoints[1] - PMTPlanePoints[0];
273  PMT_r2 = PMTPlanePoints[2] - PMTPlanePoints[0];
274  PMT_norm = PMT_r1.Cross(PMT_r2);
275  MirrorCenter.SetX(gp->fMirrorX);
276  MirrorCenter.SetY(gp->fMirrorY);
277  MirrorCenter.SetZ(gp->fMirrorZ);
278  cout << "MirrorCenter=(" << MirrorCenter.X() << "," << MirrorCenter.Y()
279  << "," << MirrorCenter.Z() << ")" << endl;
280  cout << "PMT_r1=(" << PMT_r1.X() << "," << PMT_r1.Y() << "," << PMT_r1.Z()
281  << ")" << endl;
282  cout << "PMT_r2=(" << PMT_r2.X() << "," << PMT_r2.Y() << "," << PMT_r2.Z()
283  << ")" << endl;
284  cout << "PMT_norm=(" << PMT_norm.X() << "," << PMT_norm.Y() << ","
285  << PMT_norm.Z() << ")" << endl;
286  }
287  // //HitsAndPointsWithRef();
288  HitsAndPoints();
289  RingParameters();
290  FillMcHist();
291  }
292  //fGP.Print();
293  // cout<<" ========== PMT_TRAZ = "<<TString(gSystem->Getenv("PMT_TRAZ")).Atof()<<endl;
294  // cout<<" ========== PMT_ROTX = "<<TString(gSystem->Getenv("PMT_ROTX")).Atof()<<endl;
295 }
297 void CbmRichGeoOpt::HitsAndPoints() {
298  //***********************************************************
299  //**** points
300  //***************
301  Int_t nofPoints = fRichPoints->GetEntriesFast();
302  if (nofPoints < 0 || nofPoints > 2000) { return; }
303  H_NofPhotonsPerEv->Fill(nofPoints);
304  for (Int_t ip = 0; ip < nofPoints; ip++) {
305  TVector3 PosAtDetIn;
306  TVector3 PosAtDetOut;
307  CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(ip);
308  if (NULL == point) continue;
309  //int trackId = point->GetTrackID(); if(trackId<0) continue;
310  PosAtDetIn.SetX(point->GetX());
311  PosAtDetIn.SetY(point->GetY());
312  PosAtDetIn.SetZ(point->GetZ());
313  bool Checked =
314  CheckPointLiesOnPlane(PosAtDetIn, PMTPlanePoints[0], PMT_norm);
315  if (!Checked) continue;
316  H_PointsIn_XY->Fill(PosAtDetIn.X(), PosAtDetIn.Y());
317  if (PosAtDetIn.X() <= PMTPlaneCenterX) {
318  H_PointsIn_XY_LeftHalf->Fill(PosAtDetIn.X(), PosAtDetIn.Y());
319  }
320  if (PosAtDetIn.X() > PMTPlaneCenterX) {
321  H_PointsIn_XY_RightHalf->Fill(PosAtDetIn.X(), PosAtDetIn.Y());
322  }
323  if (PosAtDetIn.X() <= PMTPlaneThirdX) {
324  H_PointsIn_XY_Left2Thirds->Fill(PosAtDetIn.X(), PosAtDetIn.Y());
325  }
326  if (PosAtDetIn.X() > PMTPlaneThirdX) {
327  H_PointsIn_XY_RightThird->Fill(PosAtDetIn.X(), PosAtDetIn.Y());
328  }
329 
330  CbmRichGeoManager::GetInstance().RotatePoint(&PosAtDetIn, &PosAtDetOut);
331  H_PointsOut_XY->Fill(PosAtDetOut.X(), PosAtDetOut.Y());
332 
334  TVector3 MomAtPMT;
335  MomAtPMT.SetX(point->GetPx());
336  MomAtPMT.SetY(point->GetPy());
337  MomAtPMT.SetZ(point->GetPz());
338  // float MagMomAtPMT=MomAtPMT.Mag();
339 
340  Int_t PointMCTrackId = point->GetTrackID();
341  if (PointMCTrackId < 0) continue;
342  CbmMCTrack* PointTrack =
343  static_cast<CbmMCTrack*>(fMcTracks->At(PointMCTrackId));
344  if (NULL == PointTrack) continue;
345  TVector3 PointMom;
346  PointTrack->GetMomentum(PointMom);
347 
348  Int_t PointMotherId = PointTrack->GetMotherId();
349  if (PointMotherId == -1) { continue; }
350  CbmMCTrack* motherTrack =
351  static_cast<CbmMCTrack*>(fMcTracks->At(PointMotherId));
352  int pdg = TMath::Abs(motherTrack->GetPdgCode());
353  int motherId = motherTrack->GetMotherId();
354  TVector3 ElMom;
355  Double_t ElTheta = 0.;
356  if (pdg == 11 && motherId == -1) {
357  motherTrack->GetMomentum(ElMom);
358  ElTheta = ElMom.Theta() * 180 / TMath::Pi();
359  }
360  double Alpha2 = MomAtPMT.Angle(PMT_norm); //*TMath::RadToDeg();
361  double Alpha2InDeg = Alpha2 * TMath::RadToDeg();
362  if (Alpha2InDeg > 90.) { Alpha2InDeg = 180. - Alpha2InDeg; }
363  H_Alpha->Fill(Alpha2InDeg);
364  if (PosAtDetIn.X() < 0. && PosAtDetIn.Y() > 0) {
365  if (ElTheta <= 25) {
366  H_Alpha_UpLeft_RegularTheta->Fill(Alpha2InDeg);
367  H_Alpha_XYposAtDet_RegularTheta->Fill(
368  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
369  }
370  H_Alpha_UpLeft->Fill(Alpha2InDeg);
371  H_Alpha_XYposAtDet->Fill(PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
372 
373  if (PosAtDetIn.X() <= PMTPlaneCenterX) {
374  H_Alpha_UpLeft_LeftHalf->Fill(Alpha2InDeg);
375  H_Alpha_XYposAtDet_LeftHalf->Fill(
376  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
377  }
378  if (PosAtDetIn.X() > PMTPlaneCenterX) {
379  H_Alpha_UpLeft_RightHalf->Fill(Alpha2InDeg);
380  H_Alpha_XYposAtDet_RightHalf->Fill(
381  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
382  }
383  if (PosAtDetIn.X() <= PMTPlaneThirdX) {
384  H_Alpha_UpLeft_Left2Thirds->Fill(Alpha2InDeg);
385  H_Alpha_XYposAtDet_Left2Thirds->Fill(
386  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
387  }
388  if (PosAtDetIn.X() > PMTPlaneThirdX) {
389  H_Alpha_UpLeft_RightThird->Fill(Alpha2InDeg);
390  H_Alpha_XYposAtDet_RightThird->Fill(
391  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
392  }
393  }
394 
395 
396  } //end loop points
397 
398 
399  //***********************************************************
400  //**** Hits
401  //***************
402  Int_t nofHits = fRichHits->GetEntriesFast();
403 
404  for (Int_t iH = 0; iH < nofHits; iH++) {
405  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(iH);
406  if (hit == NULL) continue;
407  Int_t pointInd = hit->GetRefId();
408  if (pointInd < 0) continue;
409  CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(pointInd);
410  if (point == NULL) continue;
411 
412  // H_NofPhotonsPerHit->Fill(hit->GetNPhotons());
413 
414  TVector3 inPos(point->GetX(), point->GetY(), point->GetZ());
415  TVector3 outPos;
416  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
417  H_Hits_XY->Fill(hit->GetX(), hit->GetY());
418 
419  if (hit->GetX() <= PMTPlaneCenterX) {
420  H_Hits_XY_LeftHalf->Fill(hit->GetX(), hit->GetY());
421  }
422  if (hit->GetX() > PMTPlaneCenterX) {
423  H_Hits_XY_RightHalf->Fill(hit->GetX(), hit->GetY());
424  }
425  if (hit->GetX() <= PMTPlaneThirdX) {
426  H_Hits_XY_Left2Thirds->Fill(hit->GetX(), hit->GetY());
427  }
428  if (hit->GetX() > PMTPlaneThirdX) {
429  H_Hits_XY_RightThird->Fill(hit->GetX(), hit->GetY());
430  }
431 
432  H_DiffXhit->Fill(hit->GetX() - outPos.X());
433  H_DiffYhit->Fill(hit->GetY() - outPos.Y());
434  }
435 }
436 
438 void CbmRichGeoOpt::HitsAndPointsWithRef() {
439  Int_t nofRefPoints = fRefPoints->GetEntriesFast();
440  Int_t nofPoints = fRichPoints->GetEntriesFast();
441  if (nofPoints == 0 || nofRefPoints == 0) { return; }
442  if (nofPoints > 2000) { return; }
443 
444  for (int i = 0; i < nofRefPoints; i++) {
445  TVector3 PosAtRefl;
446  TVector3 PosAtDetIn;
447  CbmRichPoint* RefPoint = (CbmRichPoint*) fRefPoints->At(i);
448  TVector3 MomAtRef;
449  MomAtRef.SetX(RefPoint->GetPx());
450  MomAtRef.SetY(RefPoint->GetPy());
451  MomAtRef.SetZ(RefPoint->GetPz()); //RefPoint->GetMomentum(MomAtRef);
452 
453  if (RefPoint == NULL) continue;
454  int RefPointTrackId = RefPoint->GetTrackID();
455  if (RefPointTrackId < 0) { continue; }
456  CbmMCTrack* RefPointTrack =
457  static_cast<CbmMCTrack*>(fMcTracks->At(RefPointTrackId));
458  int pdg0 = RefPointTrack->GetPdgCode();
459  if (TMath::Abs(pdg0) == 11) { continue; }
460  RefPoint->Position(PosAtRefl);
461  int Zpos =
462  int(10.
463  * PosAtRefl
464  .Z()); //3037 0r 3038 -->take 3038 which is the entrance point
465  //of the REFLECTED photon into the sensitive plane
466  cout << PosAtRefl.Z() << " " << Zpos << endl;
467  //if(Zpos==3037){continue;}
468  TVector3 momRef;
469  RefPointTrack->GetMomentum(momRef);
470 
471  CbmRichPoint* point = GetPMTPoint(RefPointTrackId); //
472  PosAtDetIn.SetX(point->GetX());
473  PosAtDetIn.SetY(point->GetY());
474  PosAtDetIn.SetZ(point->GetZ());
475  TVector3 MomAtPMT;
476  MomAtPMT.SetX(point->GetPx());
477  MomAtPMT.SetY(point->GetPy());
478  MomAtPMT.SetZ(point->GetPz());
479  float MagMomAtPMT = MomAtPMT.Mag();
480  //point->GetMomentum(MomAtPMT);
481 
482  Int_t PointMCTrackId = point->GetTrackID();
483  CbmMCTrack* PointTrack =
484  static_cast<CbmMCTrack*>(fMcTracks->At(PointMCTrackId));
485  if (NULL == PointTrack) continue;
486  TVector3 PointMom;
487  PointTrack->GetMomentum(PointMom);
488 
489  Int_t PointMotherId = PointTrack->GetMotherId();
490  if (PointMotherId == -1) { continue; }
491 
492  CbmMCTrack* motherTrack =
493  static_cast<CbmMCTrack*>(fMcTracks->At(PointMotherId));
494  int pdg = TMath::Abs(motherTrack->GetPdgCode());
495  int motherId = motherTrack->GetMotherId();
496  TVector3 ElMom;
497  Double_t ElTheta = 0.; //float ElMomentum=0.;
498  if (pdg == 11 && motherId == -1) {
499  motherTrack->GetMomentum(ElMom);
500  ElTheta =
501  ElMom.Theta() * 180 / TMath::Pi(); //ElMomentum=motherTrack->GetP();
502  //cout<<"ElTheta = "<<ElTheta<<", ElMomentum = "<<ElMomentum<<endl;
503  }
504 
506  bool Checked =
507  CheckPointLiesOnPlane(PosAtDetIn, PMTPlanePoints[0], PMT_norm);
508  if (!Checked)
509  continue; //cout<<" point not on plane: ("<<point->GetX()<<","<<point->GetY()<<","<<point->GetZ()<<")"<<endl; continue;
510 
511  TVector3 LineSensToPMT = PosAtDetIn - PosAtRefl;
512  float MagLineSensToPMT = LineSensToPMT.Mag();
513  H_Diff_LineRefPMT_MomAtPMT->Fill(MagLineSensToPMT - MagMomAtPMT);
514 
515  H_Theta_TwoVectors->Fill(LineSensToPMT.Angle(MomAtPMT));
517  double Alpha = LineSensToPMT.Angle(PMT_norm); //*TMath::RadToDeg();
518  double AlphaInDeg = Alpha * TMath::RadToDeg();
519  if (AlphaInDeg > 90.) { AlphaInDeg = 180. - AlphaInDeg; }
521  //double Alpha2=PointMom.Angle(PMT_norm);//*TMath::RadToDeg();
522  double Alpha2 = MomAtPMT.Angle(PMT_norm); //*TMath::RadToDeg();
523  double Alpha2InDeg = Alpha2 * TMath::RadToDeg();
524  if (Alpha2InDeg > 90.) { Alpha2InDeg = 180. - Alpha2InDeg; }
525  //cout<<PointMom.X()<<" "<<MomAtPMT.X()<<" "<<MomAtRef.X()<<" "<<Alpha<<" "<<Alpha2<<endl;
526 
527  //PosAtDetOut
528 
529 
530  H_Alpha->Fill(Alpha2InDeg);
531  if (PosAtDetIn.X() < 0. && PosAtDetIn.Y() > 0) {
532  if (ElTheta <= 25) {
533  H_Alpha_UpLeft_RegularTheta->Fill(Alpha2InDeg);
534  H_Alpha_XYposAtDet_RegularTheta->Fill(
535  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
536  }
537  H_Alpha_UpLeft->Fill(Alpha2InDeg);
538  H_Alpha_UpLeft_II->Fill(AlphaInDeg);
539  H_Alpha_XYposAtDet->Fill(PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
540 
541  if (PosAtDetIn.X() <= PMTPlaneCenterX) {
542  H_Alpha_UpLeft_LeftHalf->Fill(Alpha2InDeg);
543  H_Alpha_XYposAtDet_LeftHalf->Fill(
544  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
545  }
546  if (PosAtDetIn.X() > PMTPlaneCenterX) {
547  H_Alpha_UpLeft_RightHalf->Fill(Alpha2InDeg);
548  H_Alpha_XYposAtDet_RightHalf->Fill(
549  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
550  }
551  if (PosAtDetIn.X() <= PMTPlaneThirdX) {
552  H_Alpha_UpLeft_Left2Thirds->Fill(Alpha2InDeg);
553  H_Alpha_XYposAtDet_Left2Thirds->Fill(
554  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
555  }
556  if (PosAtDetIn.X() > PMTPlaneThirdX) {
557  H_Alpha_UpLeft_RightThird->Fill(Alpha2InDeg);
558  H_Alpha_XYposAtDet_RightThird->Fill(
559  PosAtDetIn.X(), PosAtDetIn.Y(), Alpha2InDeg);
560  }
561  }
562  }
563 
564  //***********************************************************
565  Int_t nofHits = fRichHits->GetEntriesFast();
566 
567  for (Int_t iH = 0; iH < nofHits; iH++) {
568  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(iH);
569  if (hit == NULL) continue;
570  Int_t pointInd = hit->GetRefId();
571  if (pointInd < 0) continue;
572  CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(pointInd);
573  if (point == NULL) continue;
574 
575  // H_NofPhotonsPerHit->Fill(hit->GetNPhotons());
576 
577  TVector3 inPos(point->GetX(), point->GetY(), point->GetZ());
578  TVector3 outPos;
579  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
580  H_Hits_XY->Fill(hit->GetX(), hit->GetY());
581 
582  if (hit->GetX() <= PMTPlaneCenterX) {
583  H_Hits_XY_LeftHalf->Fill(hit->GetX(), hit->GetY());
584  }
585  if (hit->GetX() > PMTPlaneCenterX) {
586  H_Hits_XY_RightHalf->Fill(hit->GetX(), hit->GetY());
587  }
588  if (hit->GetX() <= PMTPlaneThirdX) {
589  H_Hits_XY_Left2Thirds->Fill(hit->GetX(), hit->GetY());
590  }
591  if (hit->GetX() > PMTPlaneThirdX) {
592  H_Hits_XY_RightThird->Fill(hit->GetX(), hit->GetY());
593  }
594 
595  H_DiffXhit->Fill(hit->GetX() - outPos.X());
596  H_DiffYhit->Fill(hit->GetY() - outPos.Y());
597  }
598 }
599 
601 void CbmRichGeoOpt::RingParameters() {
602  Int_t nofRings = fRichRings->GetEntriesFast();
603  for (Int_t iR = 0; iR < nofRings; iR++) {
604  CbmRichRing* ring = (CbmRichRing*) fRichRings->At(iR);
605  if (NULL == ring) continue;
606  CbmTrackMatchNew* ringMatch = (CbmTrackMatchNew*) fRichRingMatches->At(iR);
607  if (NULL == ringMatch) {
608  // H_NofRings->SetBinContent(8,H_NofRings->GetBinCenter(8)+1);
609  continue;
610  }
611 
612  Int_t mcTrackId = ringMatch->GetMatchedLink().GetIndex();
613  if (mcTrackId < 0) {
614  //H_NofRings->SetBinContent(8,H_NofRings->GetBinCenter(8)+1);
615  continue;
616  } //{ continue;}
617  CbmMCTrack* mcTrack = (CbmMCTrack*) fMcTracks->At(mcTrackId);
618  if (!mcTrack) {
619  //H_NofRings->SetBinContent(8,H_NofRings->GetBinCenter(8)+1);
620  continue;
621  } // continue;
622 
623  Int_t motherId = mcTrack->GetMotherId();
624  Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
625  if (pdg != 11 || motherId != -1) {
626  //H_NofRings->SetBinContent(8,H_NofRings->GetBinCenter(8)+1);
627  continue;
628  } // continue; // only primary electrons
629 
630  Double_t momentum = mcTrack->GetP();
631  Double_t pt = mcTrack->GetPt();
632  Double_t rapidity = mcTrack->GetRapidity();
633 
634  TVector3 mom;
635  mcTrack->GetMomentum(mom);
636  Double_t theta = mom.Theta() * 180 / TMath::Pi();
637  H_MomRing->Fill(momentum);
638  H_Mom_Theta_Rec->Fill(momentum, theta);
639  H_Pt_Theta_Rec->Fill(pt, theta);
640 
641  // Double_t theta = mcTrack->Theta();
642  // cout<<" **************************** "<<theta<<endl;
643 
644  H_NofRings->Fill(nofRings);
645  if (ring->GetNofHits() >= fMinNofHits) {
646  H_Mom_Theta_Acc->Fill(momentum, theta);
647  H_Pt_Theta_Acc->Fill(pt, theta);
648  H_acc_mom_el->Fill(momentum);
649  if (theta <= 25) { H_acc_mom_el_RegularTheta->Fill(momentum); }
650 
651  H_acc_pty_el->Fill(rapidity, pt);
652  }
653 
655  float radius = ring->GetRadius();
656  if (radius <= 0.) {
657  continue;
658  } //with ideal finder --> many rings with radius -1.
659  //Test if radius is a NAN:
660  if (!(radius <= 1. || radius > 1.)) { continue; }
661  float aA = ring->GetAaxis();
662  float bA = ring->GetBaxis();
663 
664  H_Radius->Fill(radius);
665  H_aAxis->Fill(aA);
666  H_bAxis->Fill(bA);
667  H_boa->Fill(bA / aA);
668 
669  float CentX = ring->GetCenterX();
670  float CentY = ring->GetCenterY();
671 
672  H_RingCenter->Fill(CentX, CentY);
673  H_RingCenter_Aaxis->Fill(CentX, CentY, aA);
674  H_RingCenter_Baxis->Fill(CentX, CentY, bA);
675  //cout << "PMT size in x and y [cm]: " << fGP.fPmtWidthX << " " << fGP.fPmtWidthY << endl;
676  if (theta <= 25 && theta >= 24.9) {
677  H_Mom_XY_Theta25->Fill(CentX, CentY, momentum);
678  }
679 
680  H_RingCenter_boa->Fill(CentX, CentY, bA / aA);
681  if (theta <= 25) {
682  H_boa_RegularTheta->Fill(bA / aA);
683  H_RingCenter_boa_RegularTheta->Fill(CentX, CentY, bA / aA);
684  }
685  if (CentX > PMTPlaneCenterX) {
686  H_boa_RightHalf->Fill(bA / aA);
687  H_RingCenter_boa_RightHalf->Fill(CentX, CentY, bA / aA);
688  }
689  if (CentX <= PMTPlaneCenterX) {
690  H_boa_LeftHalf->Fill(bA / aA);
691  H_RingCenter_boa_LeftHalf->Fill(CentX, CentY, bA / aA);
692  }
693  if (CentX > PMTPlaneThirdX) {
694  H_boa_RightThird->Fill(bA / aA);
695  H_RingCenter_boa_RightThird->Fill(CentX, CentY, bA / aA);
696  }
697  if (CentX <= PMTPlaneThirdX) {
698  H_boa_Left2Thirds->Fill(bA / aA);
699  H_RingCenter_boa_Left2Thirds->Fill(CentX, CentY, bA / aA);
700  }
701 
702 
703  // if(CentX <= PMTPlaneCenterX && CentY >PMTPlaneCenterY){H_boa_Q1->Fill(bA/aA);}
704  int nAllHitsInR = ring->GetNofHits();
705  H_NofHitsAll->Fill(nAllHitsInR);
706  H_NofRings_NofHits->Fill(nofRings, nAllHitsInR);
707 
708  for (int iH = 0; iH < nAllHitsInR; iH++) {
709  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(ring->GetHit(iH));
710  double xH = hit->GetX();
711  double yH = hit->GetY();
712  double dRaa = aA
713  - TMath::Sqrt((CentX - xH) * (CentX - xH)
714  + (CentY - yH) * (CentY - yH));
715  H_dR_aa->Fill(dRaa);
716  double dR = radius
717  - TMath::Sqrt((CentX - xH) * (CentX - xH)
718  + (CentY - yH) * (CentY - yH));
719  H_dR->Fill(dR);
720  H_RingCenter_dR->Fill(CentX, CentY, dR);
721  if (theta <= 25) {
722  H_dR_RegularTheta->Fill(dR);
723  H_RingCenter_dR_RegularTheta->Fill(CentX, CentY, dR);
724  }
725  if (CentX > PMTPlaneCenterX) {
726  H_dR_RightHalf->Fill(dR);
727  H_RingCenter_dR_RightHalf->Fill(CentX, CentY, dR);
728  }
729  if (CentX <= PMTPlaneCenterX) {
730  H_dR_LeftHalf->Fill(dR);
731  H_RingCenter_dR_LeftHalf->Fill(CentX, CentY, dR);
732  }
733  if (CentX > PMTPlaneThirdX) {
734  H_dR_RightThird->Fill(dR);
735  H_RingCenter_dR_RightThird->Fill(CentX, CentY, dR);
736  }
737  if (CentX <= PMTPlaneThirdX) {
738  H_dR_Left2Thirds->Fill(dR);
739  H_RingCenter_dR_Left2Thirds->Fill(CentX, CentY, dR);
740  }
741  }
742  }
743 }
744 
746 void CbmRichGeoOpt::FillMcHist() {
747  for (Int_t i = 0; i < fMcTracks->GetEntriesFast(); i++) {
748  CbmMCTrack* mcTrack = (CbmMCTrack*) fMcTracks->At(i);
749  if (!mcTrack) continue;
750  Int_t motherId = mcTrack->GetMotherId();
751  Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
752 
753  TVector3 mom;
754  mcTrack->GetMomentum(mom);
755  Double_t theta = mom.Theta() * 180 / TMath::Pi();
756 
757  if (pdg == 11 && motherId == -1) {
758 
759  H_MomPt->Fill(mcTrack->GetP(), mcTrack->GetPt());
760  H_MomPrim->Fill(mcTrack->GetP());
761  H_PtPrim->Fill(mcTrack->GetPt());
762  H_Mom_Theta_MC->Fill(mcTrack->GetP(), theta);
763  H_Pt_Theta_MC->Fill(mcTrack->GetPt(), theta);
764  if (theta <= 25) { H_MomPrim_RegularTheta->Fill(mcTrack->GetP()); }
765  }
766  }
767 }
768 
770 void CbmRichGeoOpt::InitHistograms() {
771  // int nBinsX = 28; double xMin = -110.; double xMax = 110.;
772  // int nBinsY = 40; double yMin = -200; double yMax = 200.;
773 
774  H_Diff_LineRefPMT_MomAtPMT =
775  new TH1D("H_Diff_LineRefPMT_MomAtPMT",
776  "H_Diff_LineRefPMT_MomAtPMT;#Delta [cm]; Yield",
777  100,
778  -10.,
779  10.);
780 
781  H_Theta_TwoVectors = new TH1D("H_Theta_TwoVectors",
782  "H_Theta_TwoVectors;#theta [deg]; Yield",
783  100,
784  0.,
785  10.);
786  H_MomRing = new TH1D("H_MomRing", "H_MomRing;p [GeV]; Yield", 49, 0., 12.);
787  H_MomPrim = new TH1D("H_MomPrim", "H_MomPrim;p [GeV]; Yield", 49, 0., 12.);
788  H_PtPrim = new TH1D("H_PtPrim", "H_PtPrim;p [GeV]; Yield", 81, 0., 4.);
789  H_MomPt = new TH2D(
790  "H_MomPt", "H_MomPt;p [GeV];pt [GeV]; Yield", 101, 0., 10., 81, 0., 4.);
791  H_Mom_Theta_MC = new TH2D("H_Mom_Theta_MC",
792  "H_Mom_Theta_MC;p [GeV];theta [deg]; Yield",
793  40,
794  0.,
795  10.,
796  50,
797  0,
798  25.);
799  H_Pt_Theta_MC = new TH2D("H_Pt_Theta_MC",
800  "H_Pt_Theta_MC;pt [GeV];theta [deg]; Yield",
801  16,
802  0.,
803  4.,
804  50,
805  0,
806  25.);
807 
808  H_Mom_Theta_Rec = new TH2D("H_Mom_Theta_Rec",
809  "H_Mom_Theta_Rec;p [GeV];theta [deg]; Yield",
810  40,
811  0.,
812  10.,
813  50,
814  0,
815  25.);
816  H_Pt_Theta_Rec = new TH2D("H_Pt_Theta_Rec",
817  "H_Pt_Theta_Rec;pt [GeV];theta [deg]; Yield",
818  16,
819  0.,
820  4.,
821  50,
822  0,
823  25.);
824  H_Mom_Theta_Acc = new TH2D("H_Mom_Theta_Acc",
825  "H_Mom_Theta_Acc;p [GeV];theta [deg]; Yield",
826  40,
827  0.,
828  10.,
829  50,
830  0,
831  25.);
832  H_Pt_Theta_Acc = new TH2D("H_Pt_Theta_Acc",
833  "H_Pt_Theta_Acc;pt [GeV];theta [deg]; Yield",
834  16,
835  0.,
836  4.,
837  50,
838  0,
839  25.);
840 
841  H_MomPrim_RegularTheta = new TH1D("H_MomPrim_RegularTheta",
842  "H_MomPrim_RegularTheta;p [GeV]; Yield",
843  49,
844  0.,
845  12.);
846  H_acc_mom_el_RegularTheta =
847  new TH1D("H_acc_mom_el_RegularTheta",
848  "H_acc_mom_el_RegularTheta;p [GeV/c];Yield",
849  49,
850  0.,
851  12.);
852 
853 
854  H_Hits_XY = new TH2D("H_Hits_XY",
855  "H_Hits_XY;X [cm];Y [cm];Counter",
856  151,
857  -150.,
858  0.,
859  351,
860  0.,
861  350.);
862  H_Hits_XY_LeftHalf = new TH2D("H_Hits_XY_LeftHalf",
863  "H_Hits_XY_LeftHalf;X [cm];Y [cm];Counter",
864  151,
865  -150.,
866  0.,
867  351,
868  0.,
869  350.);
870  H_Hits_XY_RightHalf = new TH2D("H_Hits_XY_RightHalf",
871  "H_Hits_XY_RightHalf;X [cm];Y [cm];Counter",
872  151,
873  -150.,
874  0.,
875  351,
876  0.,
877  350.);
878  H_Hits_XY_RightThird = new TH2D("H_Hits_XY_RightThird",
879  "H_Hits_XY_RightThird;X [cm];Y [cm];Counter",
880  151,
881  -150.,
882  0.,
883  351,
884  0.,
885  350.);
886  H_Hits_XY_Left2Thirds =
887  new TH2D("H_Hits_XY_Left2Thirds",
888  "H_Hits_XY_Left2Thirds;X [cm];Y [cm];Counter",
889  151,
890  -150.,
891  0.,
892  351,
893  0.,
894  350.);
895 
896  H_PointsIn_XY = new TH2D("H_PointsIn_XY",
897  "H_PointsIn_XY;X [cm];Y [cm];Counter",
898  151,
899  -150.,
900  0.,
901  251,
902  50.,
903  250.);
904  H_PointsIn_XY_LeftHalf =
905  new TH2D("H_PointsIn_XY_LeftHalf",
906  "H_PointsIn_XY_LeftHalf;X [cm];Y [cm];Counter",
907  151,
908  -150.,
909  0.,
910  251,
911  50.,
912  250.);
913  H_PointsIn_XY_RightHalf =
914  new TH2D("H_PointsIn_XY_RightHalf",
915  "H_PointsIn_XY_RightHalf;X [cm];Y [cm];Counter",
916  151,
917  -150.,
918  0.,
919  251,
920  50.,
921  250.);
922  H_PointsIn_XY_RightThird =
923  new TH2D("H_PointsIn_XY_RightThird",
924  "H_PointsIn_XY_RightThird;X [cm];Y [cm];Counter",
925  151,
926  -150.,
927  0.,
928  251,
929  50.,
930  250.);
931  H_PointsIn_XY_Left2Thirds =
932  new TH2D("H_PointsIn_XY_Left2Thirds",
933  "H_PointsIn_XY_Left2Thirds;X [cm];Y [cm];Counter",
934  151,
935  -150.,
936  0.,
937  251,
938  50.,
939  250.);
940 
941  H_PointsOut_XY = new TH2D("H_PointsOut_XY",
942  "H_PointsOut_XY;X [cm];Y [cm];Counter",
943  151,
944  -150.,
945  0.,
946  351,
947  0.,
948  350.);
949  //cout<<" init hist H_NofPhotonsPerEv"<<endl;
950  H_NofPhotonsPerEv =
951  new TH1D("H_NofPhotonsPerEv",
952  "H_NofPhotonsPerEv;Number of photons per event;Yield",
953  500,
954  0.,
955  1000.);
956  H_NofPhotonsPerHit =
957  new TH1D("H_NofPhotonsPerHit",
958  "H_NofPhotonsPerHit;Number of photons per hit;Yield",
959  10,
960  -0.5,
961  9.5);
962  H_NofPhotonsSmallerThan30 =
963  new TH1D("H_NofPhotonsSmallerThan30",
964  "H_NofPhotonsSmallerThan30 ;Number of photons;Yield",
965  10,
966  -0.5,
967  9.5);
968  H_DiffXhit = new TH1D(
969  "H_DiffXhit", "H_DiffXhit;Y_{point}-Y_{hit} [cm];Yield", 200, -1., 1.);
970  H_DiffYhit = new TH1D(
971  "H_DiffYhit", "H_DiffYhit;Y_{point}-Y_{hit} [cm];Yield", 200, -1., 1.);
972 
973  H_Alpha = new TH1D(
974  "H_Alpha", "H_Alpha;#alpha_{photon-PMT} [deg];Yield", 180, 0., 180.);
975  H_Alpha_UpLeft = new TH1D("H_Alpha_UpLeft",
976  "H_Alpha_UpLeft;#alpha_{photon-PMT} [deg];Yield",
977  180,
978  0.,
979  180.);
980  H_Alpha_UpLeft_II =
981  new TH1D("H_Alpha_UpLeft_II",
982  "H_Alpha_UpLeft_II;#alpha_{photon-PMT} [deg];Yield",
983  180,
984  0.,
985  180.);
986  H_Alpha_UpLeft_RegularTheta =
987  new TH1D("H_Alpha_UpLeft_RegularTheta",
988  "H_Alpha_UpLeft_RegularTheta;#alpha_{photon-PMT} [deg];Yield",
989  180,
990  0.,
991  180.);
992 
993  H_Alpha_UpLeft_LeftHalf =
994  new TH1D("H_Alpha_UpLeft_LeftHalf",
995  "H_Alpha_UpLeft_LeftHalf;#alpha_{photon-PMT} [deg];Yield",
996  180,
997  0.,
998  180.);
999  H_Alpha_UpLeft_RightHalf =
1000  new TH1D("H_Alpha_UpLeft_RightHalf",
1001  "H_Alpha_UpLeft_RightHalf;#alpha_{photon-PMT} [deg];Yield",
1002  180,
1003  0.,
1004  180.);
1005  H_Alpha_UpLeft_Left2Thirds =
1006  new TH1D("H_Alpha_UpLeft_Left2Thirds",
1007  "H_Alpha_UpLeft_Left2Thirds;#alpha_{photon-PMT} [deg];Yield",
1008  180,
1009  0.,
1010  180.);
1011  H_Alpha_UpLeft_RightThird =
1012  new TH1D("H_Alpha_UpLeft_RightThird",
1013  "H_Alpha_UpLeft_RightThird;#alpha_{photon-PMT} [deg];Yield",
1014  180,
1015  0.,
1016  180.);
1017 
1018  //cout<<" init hist H_Alpha_XYposAtDet"<<endl;
1019  H_Alpha_XYposAtDet = new TH3D(
1020  "H_Alpha_XYposAtDet",
1021  "H_Alpha_XYposAtDet; X [cm]; Y [cm];#alpha_{photon-PMT} [deg];Yield",
1022  151,
1023  -150.,
1024  0.,
1025  251,
1026  50,
1027  250,
1028  180,
1029  0.,
1030  180.);
1031  H_Alpha_XYposAtDet_RegularTheta =
1032  new TH3D("H_Alpha_XYposAtDet_RegularTheta",
1033  "H_Alpha_XYposAtDet_RegularTheta; X [cm]; Y "
1034  "[cm];#alpha_{photon-PMT} [deg];Yield",
1035  151,
1036  -150.,
1037  0.,
1038  251,
1039  50,
1040  250,
1041  180,
1042  0.,
1043  180.);
1044  H_Alpha_XYposAtDet_LeftHalf =
1045  new TH3D("H_Alpha_XYposAtDet_LeftHalf",
1046  "H_Alpha_XYposAtDet_LeftHalf; X [cm]; Y [cm];#alpha_{photon-PMT} "
1047  "[deg];Yield",
1048  151,
1049  -150.,
1050  0.,
1051  251,
1052  50,
1053  250,
1054  180,
1055  0.,
1056  180.);
1057  H_Alpha_XYposAtDet_RightHalf =
1058  new TH3D("H_Alpha_XYposAtDet_RightHalf",
1059  "H_Alpha_XYposAtDet_RightHalf; X [cm]; Y [cm];#alpha_{photon-PMT} "
1060  "[deg];Yield",
1061  151,
1062  -150.,
1063  0.,
1064  251,
1065  50,
1066  250,
1067  180,
1068  0.,
1069  180.);
1070  H_Alpha_XYposAtDet_Left2Thirds =
1071  new TH3D("H_Alpha_XYposAtDet_Left2Thirds",
1072  "H_Alpha_XYposAtDet_Left2Thirds; X [cm]; Y "
1073  "[cm];#alpha_{photon-PMT} [deg];Yield",
1074  151,
1075  -150.,
1076  0.,
1077  251,
1078  50,
1079  250,
1080  180,
1081  0.,
1082  180.);
1083  H_Alpha_XYposAtDet_RightThird =
1084  new TH3D("H_Alpha_XYposAtDet_RightThird",
1085  "H_Alpha_XYposAtDet_RightThird; X [cm]; Y "
1086  "[cm];#alpha_{photon-PMT} [deg];Yield",
1087  151,
1088  -150.,
1089  0.,
1090  251,
1091  50,
1092  250,
1093  180,
1094  0.,
1095  180.);
1096 
1098  H_dFocalPoint_Delta = new TH1D("H_dFocalPoint_Delta",
1099  "H_dFocalPoint_Delta;#Delta_{f} [mm];Yield",
1100  80,
1101  -20.,
1102  20.);
1103  H_dFocalPoint_Rho = new TH1D("H_dFocalPoint_Rho",
1104  "H_dFocalPoint_Rho;#rho_{f} [mm];Yield",
1105  150,
1106  50.,
1107  200.);
1108 
1109  //cout<<" init hist H_acc_mom_el"<<endl;
1110 
1112  // Detector acceptance efficiency vs. (pt,y) and p
1113  H_acc_mom_el =
1114  new TH1D("H_acc_mom_el", "H_acc_mom_el;p [GeV/c];Yield", 49, 0., 12.);
1115  H_acc_pty_el = new TH2D("H_acc_pty_el",
1116  "H_acc_pty_el;Rapidity;P_{t} [GeV/c];Yield",
1117  25,
1118  0.,
1119  4.,
1120  61,
1121  0.,
1122  10.);
1123  H_Mom_XY_Theta25 = new TH3D("H_Mom_XY_Theta25",
1124  "H_Mom_XY_Theta25",
1125  151,
1126  -150,
1127  0,
1128  351,
1129  0,
1130  350,
1131  121,
1132  0.,
1133  12.); //25, 0., 4., 61, 0., 10.);
1135 
1136  H_NofHitsAll = new TH1D(
1137  "H_NofHitsAll", "H_NofHitsAll;Nof hits in ring;Yield", 50, 0., 50.);
1138  H_NofRings =
1139  new TH1D("H_NofRings", "H_NofRings;Nof rings per event;Yield", 10, 0., 10.);
1140  H_NofRings_NofHits =
1141  new TH2D("H_NofRings_NofHits",
1142  "H_NofRings_NofHits;Nof rings per event, Nof hits per ring;Yield",
1143  10,
1144  0.,
1145  10.,
1146  50,
1147  0.,
1148  50.);
1149 
1151  H_Radius = new TH1D("H_Radius", "H_Radius", 401, 2., 6.);
1152  H_aAxis = new TH1D("H_aAxis", "H_aAxis", 401, 2., 10.);
1153  H_bAxis = new TH1D("H_bAxis", "H_bAxis", 401, 2., 10.);
1154  H_boa = new TH1D("H_boa", "H_boa", 51, 0.5, 1.);
1155  H_boa_RegularTheta =
1156  new TH1D("H_boa_RegularTheta", "H_boa_RegularTheta", 51, 0.5, 1.);
1157  H_boa_LeftHalf = new TH1D("H_boa_LeftHalf", "H_boa_LeftHalf", 51, 0.5, 1.);
1158  H_boa_RightHalf = new TH1D("H_boa_RightHalf", "H_boa_RightHalf", 51, 0.5, 1.);
1159  H_boa_Left2Thirds =
1160  new TH1D("H_boa_Left2Thirds", "H_boa_Left2Thirds", 51, 0.5, 1.);
1161  H_boa_RightThird =
1162  new TH1D("H_boa_RightThird", "H_boa_RightThird", 51, 0.5, 1.);
1163 
1164 
1165  H_dR_aa = new TH1D("H_dR_aa", "H_dR_aa", 61, -3., 3.);
1166  H_dR = new TH1D("H_dR", "H_dR", 61, -3., 3.);
1167  H_dR_RegularTheta =
1168  new TH1D("H_dR_RegularTheta", "H_dR_RegularTheta", 61, -3., 3.);
1169  H_dR_LeftHalf = new TH1D("H_dR_LeftHalf", "H_dR_LeftHalf", 61, -3., 3.);
1170  H_dR_RightHalf = new TH1D("H_dR_RightHalf", "H_dR_RightHalf", 61, -3., 3.);
1171  H_dR_Left2Thirds =
1172  new TH1D("H_dR_Left2Thirds", "H_dR_Left2Thirds", 61, -3., 3.);
1173  H_dR_RightThird = new TH1D("H_dR_RightThird", "H_dR_RightThird", 61, -3., 3.);
1174 
1175  //cout<<" init hist H_RingCenter"<<endl;
1176 
1177  H_RingCenter =
1178  new TH2D("H_RingCenter", "H_RingCenter", 201, -100., 0., 351, 0., 350.);
1179 
1180  H_RingCenter_Aaxis = new TH3D("H_RingCenter_Aaxis",
1181  "H_RingCenter_Aaxis",
1182  151,
1183  -150,
1184  0,
1185  351,
1186  0,
1187  350,
1188  80,
1189  2.,
1190  10.);
1191  H_RingCenter_Baxis = new TH3D("H_RingCenter_Baxis",
1192  "H_RingCenter_Baxis",
1193  151,
1194  -150,
1195  0,
1196  351,
1197  0,
1198  350,
1199  80,
1200  2.,
1201  10.);
1202  H_RingCenter_boa = new TH3D("H_RingCenter_boa",
1203  "H_RingCenter_boa",
1204  151,
1205  -150,
1206  0,
1207  351,
1208  0,
1209  350,
1210  51,
1211  0.5,
1212  1.);
1213  H_RingCenter_boa_RegularTheta = new TH3D("H_RingCenter_boa_RegularTheta",
1214  "H_RingCenter_boa_RegularTheta",
1215  151,
1216  -150,
1217  0,
1218  351,
1219  0,
1220  350,
1221  51,
1222  0.5,
1223  1.);
1224  H_RingCenter_boa_LeftHalf = new TH3D("H_RingCenter_boa_LeftHalf",
1225  "H_RingCenter_boa_LeftHalf",
1226  151,
1227  -150,
1228  0,
1229  351,
1230  0,
1231  350,
1232  51,
1233  0.5,
1234  1.);
1235  H_RingCenter_boa_RightHalf = new TH3D("H_RingCenter_boa_RightHalf",
1236  "H_RingCenter_boa_RightHalf",
1237  151,
1238  -150,
1239  0,
1240  351,
1241  0,
1242  350,
1243  51,
1244  0.5,
1245  1.);
1246  H_RingCenter_boa_Left2Thirds = new TH3D("H_RingCenter_boa_Left2Thirds",
1247  "H_RingCenter_boa_Left2Thirds",
1248  151,
1249  -150,
1250  0,
1251  351,
1252  0,
1253  350,
1254  51,
1255  0.5,
1256  1.);
1257  H_RingCenter_boa_RightThird = new TH3D("H_RingCenter_boa_RightThird",
1258  "H_RingCenter_boa_RightThird",
1259  151,
1260  -150,
1261  0,
1262  351,
1263  0,
1264  350,
1265  51,
1266  0.5,
1267  1.);
1268 
1269  H_RingCenter_dR = new TH3D("H_RingCenter_dR",
1270  "H_RingCenter_dR",
1271  151,
1272  -150,
1273  0,
1274  351,
1275  0,
1276  350,
1277  61,
1278  -3.,
1279  3.);
1280  H_RingCenter_dR_RegularTheta = new TH3D("H_RingCenter_dR_RegularTheta",
1281  "H_RingCenter_dR_RegularTheta",
1282  151,
1283  -150,
1284  0,
1285  351,
1286  0,
1287  350,
1288  61,
1289  -3.,
1290  3.);
1291  H_RingCenter_dR_LeftHalf = new TH3D("H_RingCenter_dR_LeftHalf",
1292  "H_RingCenter_dR_LeftHalf",
1293  151,
1294  -150,
1295  0,
1296  351,
1297  0,
1298  350,
1299  61,
1300  -3.,
1301  3.);
1302  H_RingCenter_dR_RightHalf = new TH3D("H_RingCenter_dR_RightHalf",
1303  "H_RingCenter_dR_RightHalf",
1304  151,
1305  -150,
1306  0,
1307  351,
1308  0,
1309  350,
1310  61,
1311  -3.,
1312  3.);
1313  H_RingCenter_dR_Left2Thirds = new TH3D("H_RingCenter_dR_Left2Thirds",
1314  "H_RingCenter_dR_Left2Thirds",
1315  151,
1316  -150,
1317  0,
1318  351,
1319  0,
1320  350,
1321  61,
1322  -3.,
1323  3.);
1324  H_RingCenter_dR_RightThird = new TH3D("H_RingCenter_dR_RightThird",
1325  "H_RingCenter_dR_RightThird",
1326  151,
1327  -150,
1328  0,
1329  351,
1330  0,
1331  350,
1332  61,
1333  -3.,
1334  3.);
1335 }
1337 void CbmRichGeoOpt::WriteHistograms() {
1338 
1339  H_Diff_LineRefPMT_MomAtPMT->Write();
1340  H_Theta_TwoVectors->Write();
1341  H_MomRing->Write();
1342  H_MomPrim->Write();
1343  H_PtPrim->Write();
1344  H_MomPt->Write();
1345 
1346  H_Mom_Theta_MC->Write();
1347  H_Pt_Theta_MC->Write();
1348  H_Mom_Theta_Rec->Write();
1349  H_Pt_Theta_Rec->Write();
1350  H_Mom_Theta_Acc->Write();
1351  H_Pt_Theta_Acc->Write();
1352 
1353  H_Mom_XY_Theta25->Write();
1354 
1355  H_MomPrim_RegularTheta->Write();
1356  H_acc_mom_el_RegularTheta->Write();
1357  H_Hits_XY->Write();
1358  H_Hits_XY_LeftHalf->Write();
1359  H_Hits_XY_RightHalf->Write();
1360  H_Hits_XY_RightThird->Write();
1361  H_Hits_XY_Left2Thirds->Write();
1362 
1363  H_PointsIn_XY->Write();
1364  H_PointsIn_XY_LeftHalf->Write();
1365  H_PointsIn_XY_RightHalf->Write();
1366  H_PointsIn_XY_RightThird->Write();
1367  H_PointsIn_XY_Left2Thirds->Write();
1368 
1369  H_PointsOut_XY->Write();
1370  H_NofPhotonsPerEv->Write();
1371  H_NofPhotonsPerHit->Write();
1372  H_NofPhotonsSmallerThan30->Write();
1373  H_DiffXhit->Write();
1374  H_DiffYhit->Write();
1375 
1376  H_Alpha->Write();
1377  H_Alpha_UpLeft->Write();
1378  H_Alpha_UpLeft_II->Write();
1379  H_Alpha_UpLeft_RegularTheta->Write();
1380  H_Alpha_UpLeft_LeftHalf->Write();
1381  H_Alpha_UpLeft_RightHalf->Write();
1382  H_Alpha_UpLeft_Left2Thirds->Write();
1383  H_Alpha_UpLeft_RightThird->Write();
1384 
1385  H_Alpha_XYposAtDet->Write();
1386  H_Alpha_XYposAtDet_RegularTheta->Write();
1387  H_Alpha_XYposAtDet_LeftHalf->Write();
1388  H_Alpha_XYposAtDet_RightHalf->Write();
1389  H_Alpha_XYposAtDet_Left2Thirds->Write();
1390  H_Alpha_XYposAtDet_RightThird->Write();
1391 
1392 
1393  H_acc_mom_el->Write();
1394  H_acc_pty_el->Write();
1395 
1396  H_dFocalPoint_Delta->Write();
1397  H_dFocalPoint_Rho->Write();
1398  H_NofHitsAll->Write();
1399  H_NofRings->Write();
1400  H_NofRings_NofHits->Write();
1401  H_Radius->Write();
1402  H_aAxis->Write();
1403  H_bAxis->Write();
1404 
1405  H_boa->Write();
1406  H_boa_RegularTheta->Write();
1407  H_boa_LeftHalf->Write();
1408  H_boa_RightHalf->Write();
1409  H_boa_Left2Thirds->Write();
1410  H_boa_RightThird->Write();
1411 
1412  H_dR_aa->Write();
1413  H_dR->Write();
1414  H_dR_RegularTheta->Write();
1415  H_dR_LeftHalf->Write();
1416  H_dR_RightHalf->Write();
1417  H_dR_Left2Thirds->Write();
1418  H_dR_RightThird->Write();
1419 
1420  H_RingCenter->Write();
1421  H_RingCenter_Aaxis->Write();
1422  H_RingCenter_Baxis->Write();
1423 
1424  H_RingCenter_boa->Write();
1425  H_RingCenter_boa_RegularTheta->Write();
1426  H_RingCenter_boa_LeftHalf->Write();
1427  H_RingCenter_boa_RightHalf->Write();
1428  H_RingCenter_boa_Left2Thirds->Write();
1429  H_RingCenter_boa_RightThird->Write();
1430 
1431  H_RingCenter_dR->Write();
1432  H_RingCenter_dR_RegularTheta->Write();
1433  H_RingCenter_dR_LeftHalf->Write();
1434  H_RingCenter_dR_RightHalf->Write();
1435  H_RingCenter_dR_Left2Thirds->Write();
1436  H_RingCenter_dR_RightThird->Write();
1437 }
1440 CbmRichPoint* CbmRichGeoOpt::GetPMTPoint(int TrackIdOfSensPlane) {
1441  Int_t nofPoints = fRichPoints->GetEntriesFast();
1442  for (Int_t ip = 0; ip < nofPoints; ip++) {
1443  CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(ip);
1444  if (NULL == point) continue;
1445  int trackId = point->GetTrackID();
1446  if (trackId < 0) continue;
1447  if (
1448  trackId
1449  == TrackIdOfSensPlane) { //cout<<"In Function: got corresponding trackid:"<<trackId<<endl;
1450  return point;
1451  }
1452  }
1453  return NULL;
1454 }
1456 void CbmRichGeoOpt::FillPointsAtPMT() {
1457 
1458  for (unsigned int p = 0; p < PMTPlanePoints.size(); p++) {
1459  if (PMTPlanePoints[p].X() != -1000.) {
1460  if (p == 0) {
1461  continue;
1462  } else {
1463  int PointFilled = 1;
1464  for (int p2 = p - 1; p2 > -1; p2--) {
1465  if (TMath::Abs(PMTPlanePoints[p2].X() - PMTPlanePoints[p].X())
1466  < 1.0) {
1467  PointFilled = 0;
1468  }
1469  }
1470  if (PointFilled == 1) { continue; }
1471  }
1472  }
1473  Int_t nofPoints = fRichPoints->GetEntriesFast();
1474  // cout<<"We have "<<nofPoints<<" points registered on PMT plane"<<endl;
1475  for (Int_t ip = 0; ip < nofPoints - 10; ip += 10) {
1476  CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(ip);
1477  if (NULL == point) continue;
1478  int trackId = point->GetTrackID();
1479  if (trackId < 0) continue;
1480  if (point->GetX() >= 0 || point->GetY() <= 0) { continue; }
1481  cout << "FillPointsAtPMT: Fillinf point #" << p << " --> "
1482  << point->GetX() << ", " << point->GetY() << ", " << point->GetZ()
1483  << endl;
1484  PMTPlanePoints[p].SetX(point->GetX());
1485  PMTPlanePoints[p].SetY(point->GetY());
1486  PMTPlanePoints[p].SetZ(point->GetZ());
1487  if (PMTPlanePoints[p].X() != -1000.) { break; }
1488  }
1489  }
1490 }
1492 void CbmRichGeoOpt::FillPointsAtPMTSensPlane() {
1493 
1494  for (unsigned int p = 0; p < SensPlanePoints.size(); p++) {
1495  if (SensPlanePoints[p].X() != -1000.) {
1496  if (p == 0) {
1497  continue;
1498  } else {
1499  int PointFilled = 1;
1500  for (int p2 = p - 1; p2 > -1; p2--) {
1501  if (TMath::Abs(SensPlanePoints[p2].X() - SensPlanePoints[p].X())
1502  < 1.0) {
1503  PointFilled = 0;
1504  }
1505  }
1506  if (PointFilled == 1) { continue; }
1507  }
1508  }
1509 
1510  Int_t nofRefPoints = fRefPoints->GetEntriesFast();
1511  // cout<<"We have "<<nofRefPoints<<" points registered on sens plane"<<endl;
1512  for (Int_t ip = 0; ip < nofRefPoints - 10; ip += 10) {
1513  CbmRichPoint* RefPoint = (CbmRichPoint*) fRefPoints->At(ip);
1514  if (RefPoint == NULL) continue;
1515  int RefPointTrackId = RefPoint->GetTrackID();
1516  if (RefPointTrackId < 0) { continue; }
1517  if (RefPoint->GetX() >= 0 || RefPoint->GetY() <= 0) { continue; }
1518 
1519  SensPlanePoints[p].SetX(RefPoint->GetX());
1520  SensPlanePoints[p].SetY(RefPoint->GetY());
1521  SensPlanePoints[p].SetZ(RefPoint->GetZ());
1522  if (SensPlanePoints[p].X() != -1000.) { break; }
1523  }
1524  }
1525 }
1526 
1528 
1529 
1531 float CbmRichGeoOpt::GetIntersectionPointsLS(TVector3 MirrCenter,
1532  TVector3 G_P1,
1533  TVector3 G_P2,
1534  float R) {
1535  float A = (G_P1 - MirrCenter) * (G_P1 - MirrCenter);
1536  float B = (G_P2 - G_P1) * (G_P2 - G_P1);
1537  float P = 2. * ((G_P1 - MirrCenter) * (G_P2 - G_P1)) / (B);
1538  float q = (A - R * R) / B;
1539 
1540  float t1 = -1. * P / 2. - TMath::Sqrt((P / 2.) * (P / 2.) - q);
1541  float t2 = -1. * P / 2. + TMath::Sqrt((P / 2.) * (P / 2.) - q);
1542  //cout<<"t1="<<t1<<", t2="<<t2<<endl;
1543  //Check if nan --> no intersection
1544  if (!(t1 == 1. || t1 > 1.)) { return -1.; }
1545  //cout<<"t1="<<t1<<", t2="<<t2<<endl;
1546 
1547  TVector3 IntersectP1;
1548  TVector3 IntersectP2;
1549  IntersectP1.SetX(G_P1.X() + t1 * (G_P2.X() - G_P1.X()));
1550  IntersectP1.SetY(G_P1.Y() + t1 * (G_P2.Y() - G_P1.Y()));
1551  IntersectP1.SetZ(G_P1.Z() + t1 * (G_P2.Z() - G_P1.Z()));
1552 
1553  IntersectP2.SetX(G_P1.X() + t2 * (G_P2.X() - G_P1.X()));
1554  IntersectP2.SetY(G_P1.Y() + t2 * (G_P2.Y() - G_P1.Y()));
1555  IntersectP2.SetZ(G_P1.Z() + t2 * (G_P2.Z() - G_P1.Z()));
1556 
1557  TVector3 Line1 = IntersectP1 - G_P1;
1558  float Length1 = TMath::Sqrt(Line1.X() * Line1.X() + Line1.Y() * Line1.Y()
1559  + Line1.Z() * Line1.Z());
1560  TVector3 Line2 = IntersectP2 - G_P1;
1561  float Length2 = TMath::Sqrt(Line2.X() * Line2.X() + Line2.Y() * Line2.Y()
1562  + Line2.Z() * Line2.Z());
1563 
1564  //return Length1<Length2 ? Length1 : Length2;
1565  if (Length1 < Length2) {
1566  return Length1;
1567  } else {
1568  return Length2;
1569  }
1570 }
1572 float CbmRichGeoOpt::GetDistanceMirrorCenterToPMTPoint(TVector3 PMTpoint) {
1573  float XTerm =
1574  (PMTpoint.X() - MirrorCenter.X()) * (PMTpoint.X() - MirrorCenter.X());
1575  float YTerm =
1576  (PMTpoint.Y() - MirrorCenter.Y()) * (PMTpoint.Y() - MirrorCenter.Y());
1577  float ZTerm =
1578  (PMTpoint.Z() - MirrorCenter.Z()) * (PMTpoint.Z() - MirrorCenter.Z());
1579  return TMath::Sqrt(XTerm + YTerm + ZTerm);
1580 }
1582 bool CbmRichGeoOpt::CheckPointLiesOnPlane(TVector3 Point,
1583  TVector3 p0,
1584  TVector3 norm) {
1585  double TolaratedDiff = 0.001;
1586  double ProdP0WithNorm =
1587  p0.Dot(norm); //cout<<"ProdP0WithNorm = "<<ProdP0WithNorm;
1588  double ProdPWithNorm =
1589  Point.Dot(norm); //cout<<" ProdPWithNorm = "<<ProdPWithNorm<<endl;
1590  return TMath::Abs(ProdP0WithNorm - ProdPWithNorm)
1591  <= ((TMath::Abs(ProdP0WithNorm) < TMath::Abs(ProdPWithNorm)
1592  ? TMath::Abs(ProdPWithNorm)
1593  : TMath::Abs(ProdP0WithNorm))
1594  * TolaratedDiff);
1595 }
1597 void CbmRichGeoOpt::GetPMTRotAngels() {
1598  CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
1599  RotX = gp->fPmt.fTheta;
1600  RotY = gp->fPmt.fPhi;
1601 }
1602 
1604 //void CbmRichGeoOpt::GetPlaneCenter(float rotMir, float rotX, float rotY)
1605 //{
1606 // PMTPlaneCenterX=MinX+(MaxX-MinX)/2.; PMTPlaneCenterY=MinY+(MaxY-MinY)/2.;
1607 //}
1608 
1610 bool CbmRichGeoOpt::CheckPointLiesOnSphere(TVector3 /*Point*/) { return true; }
1611 
1613 bool CbmRichGeoOpt::CheckLineIntersectsPlane(TVector3 /*Point*/) {
1614  return true;
1615 }
1617 bool CbmRichGeoOpt::CheckLineIntersectsSphere(TVector3 /*Point*/) {
1618  return true;
1619 }
1621 void CbmRichGeoOpt::Finish() {
1622  // cout<<nPhotonsNotOnPlane<<" out of "<<nTotalPhorons<<" are not on the plane("<<float(nPhotonsNotOnPlane)/float(nTotalPhorons)<<")"<<endl;
1623  // cout<<nPhotonsNotOnSphere<<" out of "<<nTotalPhorons<<" are not on the ideal sphere("<<float(nPhotonsNotOnSphere)/float(nTotalPhorons)<<")"<<endl;
1624  WriteHistograms();
1625 }
1626 
1627 
1628 ClassImp(CbmRichGeoOpt)
CbmRichGeoOpt::H_Alpha_UpLeft_RegularTheta
TH1D * H_Alpha_UpLeft_RegularTheta
Definition: CbmRichGeoOpt.h:268
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmRichGeoOpt::FillPointsAtPMTSensPlane
void FillPointsAtPMTSensPlane()
get senspoint coordinates.
Definition: CbmRichGeoOpt.cxx:1492
CbmMCTrack::GetMomentum
void GetMomentum(TVector3 &momentum) const
Definition: CbmMCTrack.h:172
CbmRichGeoOpt::PMT_r2
TVector3 PMT_r2
Definition: CbmRichGeoOpt.h:104
CbmRichGeoOpt::H_RingCenter_Baxis
TH3D * H_RingCenter_Baxis
Definition: CbmRichGeoOpt.h:312
CbmRichGeoOpt::H_dR_RightHalf
TH1D * H_dR_RightHalf
Definition: CbmRichGeoOpt.h:307
H_boa
TH1F * H_boa
Definition: Analyze_matching.h:10
CbmRichGeoOpt::H_dR_aa
TH1D * H_dR_aa
Definition: CbmRichGeoOpt.h:304
CbmRichRecGeoPar::fMirrorZ
Double_t fMirrorZ
Definition: CbmRichRecGeoPar.h:232
CbmRichGeoOpt::H_RingCenter_boa
TH3D * H_RingCenter_boa
Definition: CbmRichGeoOpt.h:314
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmRichGeoOpt::CheckPointLiesOnSphere
bool CheckPointLiesOnSphere(TVector3 Point)
Definition: CbmRichGeoOpt.cxx:1610
CbmRichGeoOpt.h
Optimization of the RICH geometry.
CbmRichRecGeoParPmt::fPhi
Double_t fPhi
Definition: CbmRichRecGeoPar.h:58
CbmRichGeoManager::GetInstance
static CbmRichGeoManager & GetInstance()
Definition: CbmRichGeoManager.h:29
CbmRichGeoOpt::H_Mom_Theta_Rec
TH2D * H_Mom_Theta_Rec
Definition: CbmRichGeoOpt.h:226
CbmRichGeoOpt::H_PointsOut_XY
TH2D * H_PointsOut_XY
Definition: CbmRichGeoOpt.h:255
CbmRichGeoOpt::H_bAxis
TH1D * H_bAxis
Definition: CbmRichGeoOpt.h:293
CbmRichGeoOpt::H_dFocalPoint_Rho
TH1D * H_dFocalPoint_Rho
Definition: CbmRichGeoOpt.h:263
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmRichGeoOpt::H_PointsIn_XY
TH2D * H_PointsIn_XY
Definition: CbmRichGeoOpt.h:248
CbmRichGeoOpt::H_DiffXhit
TH1D * H_DiffXhit
Definition: CbmRichGeoOpt.h:259
CbmRichGeoOpt::H_RingCenter_Aaxis
TH3D * H_RingCenter_Aaxis
Definition: CbmRichGeoOpt.h:311
CbmRichRecGeoPar::fPmt
CbmRichRecGeoParPmt fPmt
Definition: CbmRichRecGeoPar.h:218
CbmRichRecGeoPar::fMirrorY
Double_t fMirrorY
Definition: CbmRichRecGeoPar.h:231
CbmRichGeoOpt::H_Alpha_XYposAtDet_RegularTheta
TH3D * H_Alpha_XYposAtDet_RegularTheta
Definition: CbmRichGeoOpt.h:275
CbmRichGeoOpt::H_NofHitsAll
TH1D * H_NofHitsAll
Definition: CbmRichGeoOpt.h:283
CbmRichGeoOpt::H_Mom_Theta_Acc
TH2D * H_Mom_Theta_Acc
Definition: CbmRichGeoOpt.h:228
CbmRichGeoOpt::H_Alpha_UpLeft_II
TH1D * H_Alpha_UpLeft_II
Definition: CbmRichGeoOpt.h:267
CbmRichGeoOpt::FillMcHist
void FillMcHist()
get MC Histos (P & Pt).
Definition: CbmRichGeoOpt.cxx:746
CbmRichGeoOpt::H_acc_mom_el_RegularTheta
TH1D * H_acc_mom_el_RegularTheta
Definition: CbmRichGeoOpt.h:238
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
CbmRichGeoOpt::GetDistanceMirrorCenterToPMTPoint
float GetDistanceMirrorCenterToPMTPoint(TVector3 PMTpoint)
calculate distance between mirror center and pmt-point.
Definition: CbmRichGeoOpt.cxx:1572
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichRecGeoPar::fMirrorX
Double_t fMirrorX
Definition: CbmRichRecGeoPar.h:230
CbmRichGeoOpt::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmRichGeoOpt.cxx:168
CbmRichGeoOpt::H_Alpha_UpLeft_Left2Thirds
TH1D * H_Alpha_UpLeft_Left2Thirds
Definition: CbmRichGeoOpt.h:272
CbmRichGeoOpt::RotX
double RotX
Definition: CbmRichGeoOpt.h:101
CbmRichGeoOpt::H_Radius
TH1D * H_Radius
Definition: CbmRichGeoOpt.h:291
H_bAxis
TH1F * H_bAxis
Definition: Analyze_matching.h:9
CbmRichRing::GetNofHits
Int_t GetNofHits() const
Definition: CbmRichRing.h:40
CbmRichRing
Definition: CbmRichRing.h:17
CbmHit::GetRefId
Int_t GetRefId() const
Definition: CbmHit.h:72
CbmRichGeoOpt::H_aAxis
TH1D * H_aAxis
Definition: CbmRichGeoOpt.h:292
CbmRichGeoOpt::fMinNofHits
Int_t fMinNofHits
Definition: CbmRichGeoOpt.h:84
CbmRichGeoOpt
Optimization of the RICH geometry.
Definition: CbmRichGeoOpt.h:44
CbmRichGeoOpt::H_Alpha_XYposAtDet_Left2Thirds
TH3D * H_Alpha_XYposAtDet_Left2Thirds
Definition: CbmRichGeoOpt.h:279
CbmRichGeoOpt::H_MomPrim_RegularTheta
TH1D * H_MomPrim_RegularTheta
Definition: CbmRichGeoOpt.h:237
CbmRichRing.h
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
CbmRichGeoOpt::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: CbmRichGeoOpt.cxx:224
CbmRichGeoOpt::H_Hits_XY
TH2D * H_Hits_XY
Definition: CbmRichGeoOpt.h:241
H_Radius
TH1F * H_Radius
Definition: Analyze_matching.h:7
CbmRichGeoOpt::fMcTracks
TClonesArray * fMcTracks
Definition: CbmRichGeoOpt.h:75
CbmRichGeoOpt::CheckLineIntersectsSphere
bool CheckLineIntersectsSphere(TVector3 Point)
Check if a given line intersects a given sphere.
Definition: CbmRichGeoOpt.cxx:1617
CbmRichGeoOpt::SensPlanePoints
vector< TVector3 > SensPlanePoints
Definition: CbmRichGeoOpt.h:112
CbmRichGeoOpt::PMTPlaneCenterY
double PMTPlaneCenterY
Definition: CbmRichGeoOpt.h:107
CbmRichRecGeoParPmt::fTheta
Double_t fTheta
Definition: CbmRichRecGeoPar.h:57
CbmRichGeoManager.h
CbmRichGeoOpt::H_NofRings_NofHits
TH2D * H_NofRings_NofHits
Definition: CbmRichGeoOpt.h:285
CbmRichGeoOpt::H_Hits_XY_LeftHalf
TH2D * H_Hits_XY_LeftHalf
Definition: CbmRichGeoOpt.h:242
CbmRichGeoOpt::GetIntersectionPointsLS
float GetIntersectionPointsLS(TVector3 MirrCenter, TVector3 G_P1, TVector3 G_P2, float R)
Definition: CbmRichGeoOpt.cxx:1531
CbmRichGeoOpt::WriteHistograms
void WriteHistograms()
write histograms to a root-file.
Definition: CbmRichGeoOpt.cxx:1337
CbmRichRecGeoParPmt::fY
Double_t fY
Definition: CbmRichRecGeoPar.h:62
CbmRichRing::GetHit
UInt_t GetHit(Int_t i) const
Definition: CbmRichRing.h:42
CbmRichGeoOpt::H_RingCenter_dR_LeftHalf
TH3D * H_RingCenter_dR_LeftHalf
Definition: CbmRichGeoOpt.h:323
CbmRichGeoOpt::H_boa
TH1D * H_boa
Definition: CbmRichGeoOpt.h:295
CbmRichGeoOpt::InitHistograms
void InitHistograms()
Initialize histograms.
Definition: CbmRichGeoOpt.cxx:770
CbmRichGeoOpt::PMTPointsFilled
Int_t PMTPointsFilled
Definition: CbmRichGeoOpt.h:82
CbmRichGeoOpt::PMTPlaneThirdX
double PMTPlaneThirdX
Definition: CbmRichGeoOpt.h:108
CbmRichGeoOpt::H_boa_Left2Thirds
TH1D * H_boa_Left2Thirds
Definition: CbmRichGeoOpt.h:300
CbmRichGeoOpt::H_Hits_XY_RightThird
TH2D * H_Hits_XY_RightThird
Definition: CbmRichGeoOpt.h:246
CbmRichGeoOpt::H_acc_mom_el
TH1D * H_acc_mom_el
Definition: CbmRichGeoOpt.h:232
CbmRichGeoOpt::GetPMTPoint
CbmRichPoint * GetPMTPoint(int TrackIdOfSensPlane)
get a point coressponding to point at sens-plane.
Definition: CbmRichGeoOpt.cxx:1440
CbmRichGeoOpt::CbmRichGeoOpt
CbmRichGeoOpt()
Standard constructor.
Definition: CbmRichGeoOpt.cxx:32
CbmRichGeoOpt::H_boa_LeftHalf
TH1D * H_boa_LeftHalf
Definition: CbmRichGeoOpt.h:297
CbmRichGeoOpt::H_RingCenter
TH2D * H_RingCenter
Definition: CbmRichGeoOpt.h:290
CbmRichRecGeoPar
PMT parameters for the RICH geometry.
Definition: CbmRichRecGeoPar.h:87
CbmRichGeoOpt::H_Alpha_XYposAtDet_RightThird
TH3D * H_Alpha_XYposAtDet_RightThird
Definition: CbmRichGeoOpt.h:278
CbmRichGeoOpt::H_dR_RightThird
TH1D * H_dR_RightThird
Definition: CbmRichGeoOpt.h:308
CbmRichGeoOpt::H_RingCenter_boa_RightHalf
TH3D * H_RingCenter_boa_RightHalf
Definition: CbmRichGeoOpt.h:317
CbmRichGeoOpt::H_RingCenter_boa_LeftHalf
TH3D * H_RingCenter_boa_LeftHalf
Definition: CbmRichGeoOpt.h:316
CbmRichGeoOpt::H_RingCenter_dR
TH3D * H_RingCenter_dR
Definition: CbmRichGeoOpt.h:321
CbmRichRingLight.h
CbmRichGeoOpt::MirrorCenter
TVector3 MirrorCenter
Definition: CbmRichGeoOpt.h:99
CbmRichGeoOpt::H_RingCenter_dR_RightThird
TH3D * H_RingCenter_dR_RightThird
Definition: CbmRichGeoOpt.h:325
CbmTrackMatchNew.h
CbmRichGeoOpt::H_dR_RegularTheta
TH1D * H_dR_RegularTheta
Definition: CbmRichGeoOpt.h:305
Point
Definition: CbmGlobalTrackingTof.cxx:90
CbmRichGeoOpt::H_boa_RightThird
TH1D * H_boa_RightThird
Definition: CbmRichGeoOpt.h:299
CbmRichGeoOpt::PMTPlaneCenterX
double PMTPlaneCenterX
Definition: CbmRichGeoOpt.h:106
CbmRichGeoOpt::PMTPlaneCenter
TVector3 PMTPlaneCenter
Definition: CbmRichGeoOpt.h:96
CbmRichGeoOpt::PMTPlanePoints
vector< TVector3 > PMTPlanePoints
Definition: CbmRichGeoOpt.h:95
CbmRichGeoOpt::H_RingCenter_boa_RegularTheta
TH3D * H_RingCenter_boa_RegularTheta
Definition: CbmRichGeoOpt.h:315
CbmRichGeoOpt::CheckLineIntersectsPlane
bool CheckLineIntersectsPlane(TVector3 Point)
Check if a given line intersects a given plane.
Definition: CbmRichGeoOpt.cxx:1613
CbmRichGeoOpt::H_Alpha_UpLeft_LeftHalf
TH1D * H_Alpha_UpLeft_LeftHalf
Definition: CbmRichGeoOpt.h:269
CbmRichGeoManager::RotatePoint
void RotatePoint(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
Definition: CbmRichGeoManager.cxx:407
CbmRichGeoOpt::fRichHits
TClonesArray * fRichHits
Definition: CbmRichGeoOpt.h:77
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmRichGeoOpt::H_Mom_Theta_MC
TH2D * H_Mom_Theta_MC
Definition: CbmRichGeoOpt.h:224
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmRichGeoOpt::HitsAndPoints
void HitsAndPoints()
Definition: CbmRichGeoOpt.cxx:297
CbmRichGeoOpt::H_Hits_XY_Left2Thirds
TH2D * H_Hits_XY_Left2Thirds
Definition: CbmRichGeoOpt.h:245
CbmRichGeoOpt::H_Alpha_UpLeft
TH1D * H_Alpha_UpLeft
Definition: CbmRichGeoOpt.h:266
CbmRichGeoOpt::H_Alpha_XYposAtDet
TH3D * H_Alpha_XYposAtDet
Definition: CbmRichGeoOpt.h:274
CbmRichRing::GetAaxis
Double_t GetAaxis() const
Definition: CbmRichRing.h:83
H_aAxis
TH1F * H_aAxis
Definition: Analyze_matching.h:8
CbmRichGeoOpt::H_dFocalPoint_Delta
TH1D * H_dFocalPoint_Delta
Definition: CbmRichGeoOpt.h:262
CbmRichGeoOpt::RotY
double RotY
Definition: CbmRichGeoOpt.h:102
CbmRichGeoOpt::CheckPointLiesOnPlane
bool CheckPointLiesOnPlane(TVector3 Point, TVector3 p0, TVector3 n)
Check if a given point lies on agiven plane.
Definition: CbmRichGeoOpt.cxx:1582
CbmRichGeoOpt::H_RingCenter_boa_Left2Thirds
TH3D * H_RingCenter_boa_Left2Thirds
Definition: CbmRichGeoOpt.h:319
CbmRichGeoOpt::H_Hits_XY_RightHalf
TH2D * H_Hits_XY_RightHalf
Definition: CbmRichGeoOpt.h:243
CbmMCTrack::GetPt
Double_t GetPt() const
Definition: CbmMCTrack.h:99
CbmRichGeoOpt::H_PointsIn_XY_RightHalf
TH2D * H_PointsIn_XY_RightHalf
Definition: CbmRichGeoOpt.h:250
CbmRichGeoOpt::H_RingCenter_dR_RightHalf
TH3D * H_RingCenter_dR_RightHalf
Definition: CbmRichGeoOpt.h:324
CbmRichGeoOpt::H_PointsIn_XY_RightThird
TH2D * H_PointsIn_XY_RightThird
Definition: CbmRichGeoOpt.h:252
CbmRichGeoOpt::H_PtPrim
TH1D * H_PtPrim
Definition: CbmRichGeoOpt.h:222
CbmRichGeoOpt::H_boa_RightHalf
TH1D * H_boa_RightHalf
Definition: CbmRichGeoOpt.h:298
CbmRichGeoOpt::H_Pt_Theta_Rec
TH2D * H_Pt_Theta_Rec
Definition: CbmRichGeoOpt.h:227
CbmRichGeoOpt::H_RingCenter_dR_Left2Thirds
TH3D * H_RingCenter_dR_Left2Thirds
Definition: CbmRichGeoOpt.h:326
CbmRichGeoOpt::PMT_norm
TVector3 PMT_norm
Definition: CbmRichGeoOpt.h:105
CbmRichGeoOpt::FillPointsAtPMT
void FillPointsAtPMT()
get point coordinates.
Definition: CbmRichGeoOpt.cxx:1456
CbmRichRecGeoParPmt::fX
Double_t fX
Definition: CbmRichRecGeoPar.h:61
CbmRichGeoOpt::H_MomPrim
TH1D * H_MomPrim
Definition: CbmRichGeoOpt.h:221
CbmRichGeoOpt::H_PointsIn_XY_Left2Thirds
TH2D * H_PointsIn_XY_Left2Thirds
Definition: CbmRichGeoOpt.h:251
CbmRichRing::GetBaxis
Double_t GetBaxis() const
Definition: CbmRichRing.h:84
CbmMCTrack::GetRapidity
Double_t GetRapidity() const
Definition: CbmMCTrack.cxx:177
CbmRichRecGeoParPmt::fZ
Double_t fZ
Definition: CbmRichRecGeoPar.h:63
CbmMCTrack.h
CbmRichGeoOpt::H_Pt_Theta_MC
TH2D * H_Pt_Theta_MC
Definition: CbmRichGeoOpt.h:225
CbmRichGeoOpt::H_dR_Left2Thirds
TH1D * H_dR_Left2Thirds
Definition: CbmRichGeoOpt.h:309
CbmRichGeoOpt::H_Alpha_UpLeft_RightHalf
TH1D * H_Alpha_UpLeft_RightHalf
Definition: CbmRichGeoOpt.h:270
CbmRichGeoOpt::H_dR
TH1D * H_dR
Definition: CbmRichGeoOpt.h:303
CbmRichGeoOpt::H_dR_LeftHalf
TH1D * H_dR_LeftHalf
Definition: CbmRichGeoOpt.h:306
CbmRichGeoOpt::H_RingCenter_dR_RegularTheta
TH3D * H_RingCenter_dR_RegularTheta
Definition: CbmRichGeoOpt.h:322
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmRichGeoOpt::HitsAndPointsWithRef
void HitsAndPointsWithRef()
Definition: CbmRichGeoOpt.cxx:438
CbmRichGeoOpt::fRefPoints
TClonesArray * fRefPoints
Definition: CbmRichGeoOpt.h:76
CbmRichGeoOpt::H_MomPt
TH2D * H_MomPt
Definition: CbmRichGeoOpt.h:223
CbmRichGeoOpt::H_Pt_Theta_Acc
TH2D * H_Pt_Theta_Acc
Definition: CbmRichGeoOpt.h:229
CbmRichGeoOpt::H_acc_pty_el
TH2D * H_acc_pty_el
Definition: CbmRichGeoOpt.h:233
CbmRichGeoOpt::PMT_r1
TVector3 PMT_r1
Definition: CbmRichGeoOpt.h:103
CbmRichRing::GetRadius
Float_t GetRadius() const
Definition: CbmRichRing.h:82
CbmRichGeoOpt::H_NofRings
TH1D * H_NofRings
Definition: CbmRichGeoOpt.h:284
CbmRichRing::GetCenterY
Float_t GetCenterY() const
Definition: CbmRichRing.h:81
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
CbmRichGeoOpt::fRichRingMatches
TClonesArray * fRichRingMatches
Definition: CbmRichGeoOpt.h:79
CbmRichGeoOpt::H_Alpha_UpLeft_RightThird
TH1D * H_Alpha_UpLeft_RightThird
Definition: CbmRichGeoOpt.h:271
CbmRichGeoOpt::H_NofPhotonsPerEv
TH1D * H_NofPhotonsPerEv
Definition: CbmRichGeoOpt.h:256
CbmRichGeoOpt::GetPMTRotAngels
void GetPMTRotAngels()
Definition: CbmRichGeoOpt.cxx:1597
CbmRichHit.h
CbmRichGeoOpt::H_Mom_XY_Theta25
TH3D * H_Mom_XY_Theta25
Definition: CbmRichGeoOpt.h:234
CbmRichGeoOpt::H_PointsIn_XY_LeftHalf
TH2D * H_PointsIn_XY_LeftHalf
Definition: CbmRichGeoOpt.h:249
CbmRichGeoOpt::fRichRings
TClonesArray * fRichRings
Definition: CbmRichGeoOpt.h:78
CbmRichGeoOpt::H_Alpha_XYposAtDet_RightHalf
TH3D * H_Alpha_XYposAtDet_RightHalf
Definition: CbmRichGeoOpt.h:277
H_dR
TH1F * H_dR
Definition: Analyze_matching.h:11
CbmRichGeoOpt::H_Theta_TwoVectors
TH1D * H_Theta_TwoVectors
Definition: CbmRichGeoOpt.h:217
CbmRichPoint
Definition: CbmRichPoint.h:24
CbmRichGeoManager::fGP
CbmRichRecGeoPar * fGP
Definition: CbmRichGeoManager.h:23
CbmRichGeoOpt::H_Diff_LineRefPMT_MomAtPMT
TH1D * H_Diff_LineRefPMT_MomAtPMT
histograms.
Definition: CbmRichGeoOpt.h:216
CbmRichGeoOpt::H_Alpha_XYposAtDet_LeftHalf
TH3D * H_Alpha_XYposAtDet_LeftHalf
Definition: CbmRichGeoOpt.h:276
CbmRichGeoOpt::~CbmRichGeoOpt
virtual ~CbmRichGeoOpt()
Standard destructor.
Definition: CbmRichGeoOpt.cxx:166
CbmRichGeoOpt::H_DiffYhit
TH1D * H_DiffYhit
Definition: CbmRichGeoOpt.h:260
CbmRichGeoOpt::H_RingCenter_boa_RightThird
TH3D * H_RingCenter_boa_RightThird
Definition: CbmRichGeoOpt.h:318
CbmRichGeoOpt::H_MomRing
TH1D * H_MomRing
Definition: CbmRichGeoOpt.h:231
CbmMCTrack::GetP
Double_t GetP() const
Definition: CbmMCTrack.h:100
CbmRichGeoOpt::H_boa_RegularTheta
TH1D * H_boa_RegularTheta
Definition: CbmRichGeoOpt.h:296
CbmRichRing::GetCenterX
Float_t GetCenterX() const
Definition: CbmRichRing.h:80
CbmRichGeoOpt::fEventNum
Int_t fEventNum
Definition: CbmRichGeoOpt.h:81
CbmRichGeoOpt::fRichPoints
TClonesArray * fRichPoints
Definition: CbmRichGeoOpt.h:74
CbmRichHit
Definition: CbmRichHit.h:19
CbmRichGeoOpt::SensPointsFilled
Int_t SensPointsFilled
Definition: CbmRichGeoOpt.h:111
CbmRichGeoOpt::H_NofPhotonsPerHit
TH1D * H_NofPhotonsPerHit
Definition: CbmRichGeoOpt.h:257
CbmRichGeoOpt::H_NofPhotonsSmallerThan30
TH1D * H_NofPhotonsSmallerThan30
Definition: CbmRichGeoOpt.h:258
CbmRichGeoOpt::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmRichGeoOpt.cxx:1621
CbmRichGeoOpt::H_Alpha
TH1D * H_Alpha
Definition: CbmRichGeoOpt.h:264
CbmRichGeoOpt::RingParameters
void RingParameters()
Loop over all rings in array and fill ring parameters histograms.
Definition: CbmRichGeoOpt.cxx:601