FairRoot/PandaRoot
PndDrcHitProducerReal.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndDrcHitProducerReal source file -----
3 // ----- Created 30/10/09 by Dipanwita Dutta -----
4 // ----- -----
5 // ----- -----
6 // -------------------------------------------------------------------------
7 #include <fstream>
8 #include <iostream>
9 #include "stdio.h"
10 
11 #include "PndGeoDrc.h"
12 #include "PndDrcHitProducerReal.h"
13 #include "FairRootManager.h"
14 #include "PndMCTrack.h"
15 #include "PndDrcBarPoint.h"
16 #include "PndDrcPDPoint.h"
17 #include "PndDrcHit.h"
18 #include "PndDrcPDHit.h"
19 #include "TVector3.h"
20 #include "TRandom.h"
21 #include "FairRunAna.h"
22 #include "FairRuntimeDb.h"
23 #include "FairBaseParSet.h"
24 #include "FairGeoVolume.h"
25 #include "TString.h"
26 #include "FairGeoTransform.h"
27 #include "FairGeoVector.h"
28 #include "FairGeoMedium.h"
29 #include "FairGeoNode.h"
30 #include "PndGeoDrcPar.h"
31 #include "TFormula.h"
32 #include "TMath.h"
33 #include "TParticlePDG.h"
34 #include "TDatabasePDG.h"
35 #include "TPDGCode.h"
36 #include "TGeoManager.h"
37 #include "TF1.h"
38 #include "TF2.h"
39 #include "TH2.h"
40 
41 using std::endl;
42 using std::cout;
43 
44 // ----- Default constructor -------------------------------------------
45 PndDrcHitProducerReal::PndDrcHitProducerReal()
46  :PndPersistencyTask("PndDrcHitProducerReal"),fisDetEff(kTRUE),fisPixel(kTRUE),fGeo(new PndGeoDrc()),fVerbose(0),fDetType(1){
47  SetParameters();
48 }
49 
50 // ----- Standard constructor with verbosity level -------------------------------------------
51 PndDrcHitProducerReal::PndDrcHitProducerReal(Int_t verbose, Int_t det_type)
52  :PndPersistencyTask("PndDrcHitProducerReal"),fisDetEff(kTRUE),fisPixel(kTRUE),fGeo(new PndGeoDrc()),fVerbose(verbose),fDetType(det_type){
53  SetParameters();
54 }
55 
56 //-------------Set Parameter------------------------------------
57 void PndDrcHitProducerReal::SetParameters(){
58  fPixelDim=0.65; //3.1 Pixel Dimension of photocathode is 6.5mm
59  nRefrac=fGeo->nEV(); //1.467; //Refractive index of SOB
60  fSigmat=0.05; //Time Resolution is 50 ps ############################
61  fCollectionEff=0.65; //Collection Efficiency
62  fPackingFraction=0.80; //Packing Efficiency
63  fRoughness = 0.001; // 10 A
64 
65  // basic DIRC parameters:
66  fpi = TMath::Pi();
67  fzup = fGeo->barBoxZUp();
68  fzdown = fGeo->barBoxZDown();
69  fradius = fGeo->radius(); //radius in the middle of the bar = 50.cm
70  fhthick = fGeo->barHalfThick(); //half thickness of the bars=1.7/2 cm
71  fpipehAngle = fGeo->PipehAngle();
72  fbbGap = fGeo->BBoxGap();
73  fbbnum = fGeo->BBoxNum();
74  fbarnum = fGeo->barNum();
75  fphi0 = (180.-2.*fpipehAngle)/fbbnum + fpipehAngle;
76  fdphi = (180.-2.*fpipehAngle)/fbbnum*2.;
77  flside = fGeo->Lside();
78  fbarwidth = fGeo->BarWidth();
79 }
80 
81 // ----- Destructor ----------------------------------------------------
82 PndDrcHitProducerReal::~PndDrcHitProducerReal(){
83  if (fGeo) delete fGeo;
84 }
85 
86 // ----- Initialization -----------------------------------------------
87 InitStatus PndDrcHitProducerReal::Init()
88 {
89  cout << " ---------- INITIALIZATION ------------" << endl;
90  nevents = 0;
91  // Get RootManager
92  FairRootManager* ioman = FairRootManager::Instance();
93  if ( ! ioman ) {
94  cout << "-E- PndDrcHitProducerReal::Init: "
95  << "RootManager not instantiated!" << endl;
96  return kFATAL;
97  }
98  // Get input array
99  fMCArray = (TClonesArray*) ioman->GetObject("MCTrack");
100  if ( ! fMCArray ) {
101  cout << "-W- PndDrcRecoLookupMap::Init: "
102  << "No MCTrack array!" << endl;
103  return kERROR;
104  }
105 
106  // Get input array
107  fBarPointArray = (TClonesArray*) ioman->GetObject("DrcBarPoint");
108  if ( ! fBarPointArray ) {
109  cout << "-W- PndDrcHitProducerReal::Init: "
110  << "No DrcBarPoint array!" << endl;
111  return kERROR;
112  }
113 
114  // Get Photon point array
115  fPDPointArray = (TClonesArray*) ioman->GetObject("DrcPDPoint");
116  if ( ! fPDPointArray ) {
117  cout << "-W- PndDrcAna::Init: "
118  << "No DrcPDPoint array!" << endl;
119  return kERROR;
120  }
121 
122  // Create and register output array
123  fHitArray = new TClonesArray("PndDrcHit");
124  ioman->Register("DrcHit","Drc",fHitArray, GetPersistency());
125 
126  // Create and register output array
127  fPDHitArray = new TClonesArray("PndDrcPDHit");
128  ioman->Register("DrcPDHit","Drc",fPDHitArray, GetPersistency());
129 
130  cout << "-I- PndDrcHitProducerReal: Intialization successfull" << endl;
131 
132  return kSUCCESS;
133 
134 }
135 // ----- Execution of Task ---------------------------------------------
136 void PndDrcHitProducerReal::Exec(Option_t*) // ion //[R.K.03/2017] unused variable(s)
137 {
138  if ( ! fHitArray ) Fatal("Exec", "No HitArray");
139  fHitArray->Delete();
140  if ( ! fPDHitArray ) Fatal("Exec", "No Photon HitArray");
141  fPDHitArray->Delete();
142  nevents++;
143  if (fVerbose > 1) printf("\n\n=====> Event No. %d\n", nevents);
144  if (fVerbose > 1) std::cout<<"=====> Event No. "<<nevents<<endl;
145 
146  ProcessBarPoint();
147  ProcessPhotonPoint();
148 }
149 // ----- Process Bar Point ---------------------------------------------
150 void PndDrcHitProducerReal::ProcessBarPoint()
151 {
152  PndDrcBarPoint* pt=NULL;
153 
154  if (fVerbose > 0) {
155  cout << " ----------------- DRC Hit Producer --------------------" << endl;
156  cout <<" Number of input MC points in the bar: "<<fBarPointArray->GetEntries()<<endl;
157  }
158 
159  // fNHits = 0;
160 
161  // Loop over PndDrcPoints
162  for(Int_t j=0; j<fBarPointArray->GetEntriesFast(); j++) {
163 
164  pt = (PndDrcBarPoint*)fBarPointArray->At(j);
165 
166  Double_t Px= pt->GetPx();
167  Double_t Py= pt->GetPy();
168  Double_t Pz= pt->GetPz();
169  Double_t P = sqrt(Px*Px + Py*Py +Pz*Pz);
170  Double_t mass = pt->GetMass();
171  Double_t energy = TMath::Sqrt(P*P + mass*mass);
172  TVector3 MomBar;
173  pt->Momentum(MomBar);
174  pt->Position(fPosHit);
175  Double_t beta;
176  if(energy != 0) {
177  beta = P/energy;
178  }
179  else {
180  beta = -1.;
181  if (fVerbose >0) cout << "Beta not calculated " << endl;
182  }
183 
184  if (pt->GetThetaC() != -1. && beta > 1/1.47){
185  fDetectorID = pt->GetBoxId()*10 + pt->GetBarId();
186 
187  // calculate the center of the bars from the detectorID
188  //Int_t s = (fDetectorID /10);// correction DD //[R.K. 01/2017] unused variable
189  //Int_t b = (fDetectorID % 10); //[R.K. 01/2017] unused variable?
190 
191  //cout<<"-I- HitProducerReal: s = "<<s<<", b = "<<b<<endl;
192 
193  //Double_t bbAngle = ( 180. - 2.*fpipehAngle - fbbGap/fradius/fpi*180.*(fbbnum/2.-1.) )/(fbbnum/2.); //[R.K. 01/2017] unused variable?
194  //Double_t bbX = fradius*bbAngle/180.*fpi; //[R.K. 01/2017] unused variable?
195 
196  //Double_t phi_curr = (90. - fphi0 - fdphi*(s-1))/180.*fpi; //[R.K. 01/2017] unused variable
197  //if(s > fbbnum/2){ phi_curr = (90. - fphi0 - fdphi*(s-1) - 2.*fpipehAngle)/180.*fpi; } //[R.K. 01/2017] unused variable
198  //cout<<"-I- HitProducerReal: phi_curr = "<< phi_curr/fpi*180.<<endl;
199  //Double_t Xs = fradius * cos(phi_curr); //[R.K. 01/2017] unused variable?
200  //Double_t Ys = fradius * sin(phi_curr); //[R.K. 01/2017] unused variable?
201  //cout<<"-I- HitProducerReal: Xs = "<<Xs<<", Ys = "<<Ys<<endl;
202  //Double_t Xb = bbX/2.*((2.*b-1.)/fbarnum - 1.)*sin(phi_curr); //[R.K. 01/2017] unused variable?
203  //Double_t Yb = -bbX/2.*((2.*b-1.)/fbarnum - 1.)*cos(phi_curr); //[R.K. 01/2017] unused variable?
204  //cout<<"-I- HitProducerReal: ((2.*b-1.)/barnum - 1.)"<< ((2.*b-1.)/barnum - 1.)<<endl;
205  //cout<<"-I- HitProducerReal: Xb = "<<Xb<<", Yb = "<<Yb<<endl;
206  //Double_t fXHit = Xs+Xb; //[R.K. 01/2017] unused variable?v
207  //Double_t fYHit = Ys+Yb; //[R.K. 01/2017] unused variable?
208  //Double_t fZHit = 0.;
209  //Double_t fZHit = fGeo->barBoxZUp() - fGeo->EVlen(); //[R.K. 01/2017] unused variable?
210 
211  //cout<<"-I- HitProducerReal: bar center = "<< fXHit << ", "<< fYHit <<endl;
212 
213  //cout << "-I- HitProducerReal: hit phi: "<< acos(fXHit/radius)/pi*180. << endl;
214 
215  // fPosHit.SetXYZ(fXHit,fYHit,fZHit);
216 
217  Double_t fDPosXHit = 0.5; //mm
218  Double_t fDPosYHit = 0.5;
219  Double_t fDPosZHit = 0.;
220  fDPosHit.SetXYZ(fDPosXHit,fDPosYHit,fDPosZHit);
221  // realistic Values are changed after considering the tracking effects by H.Kumawat, 7 March 2014,(h.kumawat@gsi.de)
222  Double_t theta = MomBar.Theta()*180.0/TMath::Pi();
223  Float_t sigma=0.002;
224  if(theta<45.){
225  sigma=0.001*(3.21659E-01 + 7.48416E-02*theta - 9.87561E-05*theta*theta - 2.47129E-06*theta*theta*theta);//prism
226  }else if(theta>=45. && theta<90.){
227  sigma=0.001*(3.21659E-01 + 7.48416E-02*theta - 9.87561E-05*theta*theta - 2.47129E-06*theta*theta*theta);//prism
228  }else if(theta>=90.){
229  sigma=0.001*( 2.34224e+01-3.52791e-01*theta+1.64046e-03*theta*theta-7.57365e-07*theta*theta*theta);
230  }
231  sigma=sigma+0.00873*exp(-MomBar.Mag()/0.4614);
232  fThetaC = gRandom->Gaus(pt->GetThetaC(),sigma); // changed by H.Kumawat on 7 March 2014 (h.kumawat@gsi.de)
233  fErrThetaC = 0.; //rad
234 
235  fRefIndex = j;
236 
237 
238  AddHit(fDetectorID,
239  fPosHit,
240  fDPosHit,
241  fThetaC,
242  fErrThetaC,
243  fRefIndex);
244 
245  }
246  }
247 }
248 //-------- Photon Detector Hit production with efficiency-------------
249 void PndDrcHitProducerReal::ProcessPhotonPoint()
250 {
251  if (fVerbose > 0) {
252  cout <<" Number of Photon MC Points in Photon Detector Plane : "<<fPDPointArray->GetEntries()<<endl;
253  }
254 
255  Double_t lambda_min,lambda_max,lambda_step;
256  Double_t DetEfficiency[800];
257  //SetFakeDetEff(lambda_min,lambda_max,lambda_step,DetEfficiency);
258  SetPhotonDetEffNew(lambda_min,lambda_max,lambda_step,DetEfficiency);
259  // SetPhotonDetEff(lambda_min,lambda_max,lambda_step,DetEfficiency);
260  // SetPhotonDetEffOld(lambda_min,lambda_max,lambda_step,DetEfficiency);
261 
262  Double_t lambda_min_tr,lambda_max_tr,lambda_step_tr, angle_step_tr;
263  Int_t lambda_points_tr;
264  Double_t TransportEfficiency[798];
265  SetPhotonTransportEff(lambda_min_tr,lambda_max_tr,lambda_step_tr,angle_step_tr, lambda_points_tr,TransportEfficiency);
266 
267  PndDrcPDPoint* Ppt=NULL;
268  PndMCTrack* tr = NULL;
269 
270  // Loop over PndDrcPDPoints
271  for(Int_t k=0; k<fPDPointArray->GetEntriesFast(); k++) {
272 
273  Ppt = (PndDrcPDPoint*)fPDPointArray->At(k);
274 
275  Int_t trID= Ppt->GetTrackID();
276  tr = (PndMCTrack*)fMCArray->At(trID);
277 
278  // production point of the photon
279  TVector3 StartVertex = tr->GetStartVertex();
280  // initial momentum of the photon
281  TVector3 PphoInit;
282  PphoInit.SetXYZ(tr->GetMomentum().X(), tr->GetMomentum().Y(), tr->GetMomentum().Z());
283  //calculate the number of bounces:
284  Int_t NbouncesX, NbouncesY;
285  Double_t angleX, angleY;
286  NumberOfBounces(StartVertex, PphoInit, &NbouncesX, &NbouncesY, &angleX, &angleY);
287  //cout<<"N bouncesX = "<<NbouncesX<<", NbouncesY = "<<NbouncesY<<", angleX = "<<angleX<<", angleY = "<<angleY<<endl;
288 
289  Double_t PPx= Ppt->GetPx();
290  Double_t PPy= Ppt->GetPy();
291  Double_t PPz= Ppt->GetPz();
292 
293  Double_t etot = sqrt(PPx*PPx + PPy*PPy +PPz*PPz);// in GeV
294  Double_t lambda=197.0*2.0*fpi/(etot*1.0E9);//wavelength of photon in nm
295  //cout<<"1. lambda = "<<lambda<<endl;
296 
297  if(fisDetEff){
298  //cout<<"det eff!!! lam_min = "<<lambda_min<<", lam_max = "<<lambda_max<<endl;
299  if (lambda >= lambda_min && lambda < lambda_max) {
300  Int_t ilambda=(Int_t)((lambda-lambda_min)/lambda_step);
301  Double_t rand = gRandom->Rndm();
302  detection = 0;
303  if (DetEfficiency[ilambda]*fCollectionEff*fPackingFraction > rand ) detection = 1;
304  }
305  }
306  if(!fisDetEff){
307  detection=1;
308  }
309  //cout<<"detection = "<<detection<<endl;
310  //####################################
311  // transport efficiency
312  // assume that detection efficiency above CAN NOT be used together with transport efficiency
313  // Maria Patsyuk
314  //cout<<"do transport eff? "<<fisTransportEff<<", lambda = "<<lambda<<endl;
315  if(fisTransportEff){
316  if (lambda >= lambda_min_tr && lambda < lambda_max_tr) {
317  Int_t ilambda=(Int_t)((lambda-lambda_min_tr)/lambda_step_tr);
318  Int_t iangleX =(Int_t)(angleX/angle_step_tr);
319  Int_t iangleY =(Int_t)(angleY/angle_step_tr);
320  //cout<<"iangleX = "<<iangleX<<", iangleY = "<<iangleY<<", ilam = "<<ilambda<<", lambda_points_tr = "<<lambda_points_tr<<endl;
321  Double_t rand = gRandom->Rndm();
322  detection = 0;
323  //cout<<"Bounce eff = "<<TransportEfficiency[ilambda]<<endl;
324  Double_t TotalTrProb = 1.;
325  //cout<<"before ref prob. angleX = "<<angleX<<", angleY = "<<angleY<<endl;
326  Double_t ReflectionProbX = TransportEfficiency[ilambda+ iangleX*lambda_points_tr];
327  Double_t ReflectionProbY = TransportEfficiency[ilambda+ iangleY*lambda_points_tr];
328  TotalTrProb = pow(ReflectionProbX, (Int_t)NbouncesX)*pow(ReflectionProbY, (Int_t)NbouncesY);
329  //cout<<"tr eff X = "<<ReflectionProbX<<", tr eff Y = "<<ReflectionProbY<<", total = "<<TotalTrProb<<endl;
330  if(TotalTrProb > rand) detection = 1;
331  }
332  }
333  if(!fisTransportEff && !fisDetEff){
334  detection=1;
335  }
336  //####################################
337 
338  Double_t xP= Ppt->GetX();
339  Double_t yP= Ppt->GetY();
340  Double_t zP= Ppt->GetZ();
341 
342  Double_t xHit, yHit, zHit;
343  Int_t pmtID;
344 
345  // fPosPDHit.SetXYZ(xP,yP,zP);
346 
347  zHit=zP;
348  pmtID=k;
349  if(fisPixel){
350  if(ftilt == 0.){ // ###########################
351  FindDrcHitPosition(xP, yP, xHit, yHit, pmtID);
352  } // ###########################
353  //#############################################
354  // in case of tilted PD plane
355  // Maria Patsyuk
356  if(ftilt != 0.){
357 
358  //Int_t s = (fDetectorID /10);// correction DD //[R.K. 01/2017] unused variable?
359 
360  TVector3 point;
361  point.SetXYZ(xP,yP,0.);
362  //cout<<"-I- HitProducerReal: pointx = "<<xP<<", pointy = "<<yP<<endl;
363  Double_t phiP = point.Phi()/fpi*180.; // degrees
364  if(phiP < 0.){phiP = 360. + point.Phi()/3.1415*180.;}
365  Double_t phi_rot = 0.;
366  //cout<<"-I- HitProducerReal: phiP = "<<phiP<<endl;
367  if(phiP > 0. && phiP < 86.4){
368  phi_rot = TMath::Floor(phiP/fdphi) *fdphi + fdphi/2.;
369  }
370  if(phiP > 93.6 && phiP < 266.4){
371  phi_rot = 90. + fpipehAngle + TMath::Floor((phiP-90.-fpipehAngle)/fdphi) *fdphi + fdphi/2.;
372  }
373  if(phiP > 273.6 && phiP < 360.){
374  phi_rot = 270. + fpipehAngle + TMath::Floor((phiP-270.-fpipehAngle)/fdphi) *fdphi + fdphi/2.;
375  }
376  if(phiP > 86.4 && phiP < 93.6){
377  phi_rot = 90.;
378  }
379  if(phiP > 266.4 && phiP < 273.6){
380  phi_rot = 270.;
381  }
382  //cout<<"-I- HitProducerReal: phi_rot = "<<phi_rot<<endl;
383  TVector3 vhit;
384  vhit.SetXYZ(xP, yP, 0.);
385  vhit.RotateZ(-phi_rot/180.*fpi); // rad
386  Double_t xP_bar = -vhit.Y();
387  Double_t yP_bar = vhit.X();
388  Double_t yHit_bar, xHit_bar;
389  FindDrcHitPositionTilt(xP_bar, yP_bar, xHit_bar, yHit_bar, pmtID);
390  TVector3 hitbar;
391  hitbar.SetXYZ(yHit_bar, -xHit_bar, 0.);
392  //cout<<"-I- HitProducerReal: BAR hitx = "<<xHit_bar<<", yhit = "<<yHit_bar<<endl;
393  //cout<<"-I- HitProducerReal: phi hit bar = "<<vhit.Phi()/3.1415*180.<<", hitbar phi = "<<hitbar.Phi()/3.1415*180.<<endl;
394  //cout<<"-I- HitProducerReal: dphi = "<<fdphi<<endl;
395  if(hitbar.Phi() < fdphi/2./180.*3.1415 && hitbar.Phi() > -fdphi/2./180.*3.1415){
396  vhit.SetXYZ(yHit_bar, -xHit_bar, 0.);
397  vhit.RotateZ(phi_rot/180.*fpi);
398  //cout<<"-I- HitProducerReal: xhit = "<<vhit.X()<<", yhit = "<<vhit.Y()<<endl;
399  xHit = TMath::Nint(vhit.X()*1000.)/1000.;
400  yHit = TMath::Nint(vhit.Y()*1000.)/1000.;
401  //cout<<"-I- HitProducerReal: hitx = "<<xHit<<", yhit = "<<yHit<<endl;
402  }
403  else {
404  continue;
405  }
406  }
407  //#############################################
408 
409  }else{
410  xHit=xP;
411  yHit=yP;
412  }
413  fPosPDHit.SetXYZ(xHit,yHit,zHit);
414 
415  Double_t fDPosXPDHit = fPixelDim/2; //mm
416  Double_t fDPosYPDHit = fPixelDim/2; //mm
417  Double_t fDPosZPDHit = 0.;
418 
419  fDPosPDHit.SetXYZ(fDPosXPDHit,fDPosYPDHit,fDPosZPDHit);
420 
421  Double_t time=Ppt->GetTime();
422  Smear(time,fSigmat);
423  fTime=time;
424 
425  fTimeThreshold = 500.;
426  fPDRefIndex = k;
427 
428  if(detection==1){
429  AddPDHit(fDetectorID,
430  fPosPDHit,
431  fDPosPDHit,
432  fTime,
433  fTimeThreshold,
434  fPDRefIndex);
435  }
436  }
437 }
438 //------ Find Nubmer of Bounces --------------------------------------
439 void PndDrcHitProducerReal::NumberOfBounces(TVector3 start, TVector3 dir, Int_t *n1, Int_t *n2, Double_t *alpha1, Double_t *alpha2){
440  // calculates the number of bounces in x and y direction and reflection angles in these directions.
441 
442  Double_t PhiRot = FindPhiRot(start.X(), start.Y());
443  //cout<<"-I- NumberOfBounces: phi rot = "<<PhiRot<<endl;
444 
445  // Photon production point in bar' coordinate system (origin at the corner of the bar):
446  TVector3 startBar;
447  startBar.SetXYZ(start.X(), start.Y(), start.Z());
448  //cout<<"-I- NumberOfBounces: start.X = "<<start.X()<<", start.Y = "<<start.Y()<<endl;
449  startBar.RotateZ(-PhiRot/180.*fpi);
450  //cout<<"-I- NumberOfBounces: startBar.X = "<<startBar.X()<<", startBar.Y = "<<startBar.Y()<<endl;
451 
452  // Photon momentum in bar' coord system:
453  TVector3 PphoB;
454  PphoB = dir.Unit();
455  PphoB.RotateZ(-PhiRot/180.*fpi);
456 
457  // Find coordinates of X0, Y0:
458  Double_t Z0, X0, Y0;
459  if(dir.Theta() < 3.1415/2.){
460  Z0 = -(fabs(fzup) + 2.*fzdown - startBar.Z());
461  }
462  if(dir.Theta() >= 3.1415/2.){
463  Z0 = -(startBar.Z() - fzup);
464  }
465  //cout<<"-I- NumberOfBounces: Z0 = "<<Z0<<", Theta = "<<PphoB.Theta()/3.1415*180.<<", tan t = "<<tan(PphoB.Theta())<<", phi = "<<PphoB.Phi()/3.1415*180.<<endl;
466  X0 = Z0*tan(PphoB.Theta())*cos(PphoB.Phi());
467  Y0 = Z0*tan(PphoB.Theta())*sin(PphoB.Phi());
468  //cout<<"-I- NumberOfBounces: X0 = "<<X0<<", Y0 = "<<Y0<<endl;
469 
470  // Find the number of bounces in each direction
471  Double_t N1, N2;
472  //frad_out = (fradius-fhthick)/cos(2.*3.1415/16./2.); // radius at corner - thickness ###
473  //flside = 2.*frad_out*sin(2.*3.1415/16./2.) - (2.*fboxthick) - (2.*fboxgap);
474  //flside = (180. - 2.*fpipehAngle - fbbGap/fradius*(fbbnum/2. - 1.)/fpi*180.)/(fbbnum/2.) * fradius/ 180.*fpi;
475  //fbarwidth = flside/fbarnum;
476  //cout<<"-I- NumberOfBounces: lside = "<<flside<<", bar width = "<<fbarwidth<<endl;
477 
478  //cout<<"fbarnum = "<<fbarnum<<endl;
479  // Find which bar in the bar box was hit:
480  if(fbarnum > 1){
481  Int_t NhitBar = (Int_t)((0.5*flside + startBar.Y())/fbarwidth)+1;
482  //cout<<"-I- NumberOfBounces: bar "<<NhitBar<<" was hit, "<<((0.5*flside + startBar.Y())/fbarwidth)+1<<endl;
483  //cout<<"-I- NumberOfBounces: start bar Y = "<<startBar.Y()<<endl;
484 
485  //cout<<"-I- NumberOfBounces: start position X = "<< startBar.X() - (fradius-fhthick)<<", Y = "<<startBar.Y() + 0.5*flside-(NhitBar-1)*fbarwidth<<endl;
486  FindOutPoint(X0, startBar.X()-(fradius-fhthick), 2.*fhthick, &N1, 0);
487  FindOutPoint(Y0, startBar.Y()+0.5*flside-(NhitBar-1)*fbarwidth, fbarwidth, &N2, 0);
488  //cout<<"-I- NumberOfBounces: N1 = "<<N1<<", N2 = "<<N2<<endl;
489  }
490 
491  if(fbarnum == 1 && flside > fbarwidth){
492  //cout<<"second case!!!"<<endl;
493  FindOutPoint(X0, startBar.X()-(fradius-fhthick), 2.*fhthick, &N1, 0);
494  FindOutPoint(Y0, startBar.Y()+0.5*fbarwidth, fbarwidth, &N2, 0);
495  }
496 
497  *n1 = (Int_t)N1;
498  *n2 = (Int_t)N2;
499 
500  // calculate the reflection angles in x and y directions:
501  TVector3 up_down;
502  up_down.SetXYZ(0.,1.,0.);
503  TVector3 left_right;
504  left_right.SetXYZ(1.,0.,0.);
505  Double_t angle1 = PphoB.Angle(left_right);
506  if(angle1 > fpi/2.){angle1 = fpi - PphoB.Angle(left_right);}
507  Double_t angle2 = PphoB.Angle(up_down);
508  if(angle2 > fpi/2.){angle2 = fpi - PphoB.Angle(up_down);}
509  *alpha1 = angle1;
510  *alpha2 = angle2;
511  //cout<<"-I- NumberOfBounces: angle1 = "<<angle1<<", angle2 = "<<angle2<<endl;
512 }
513 
514 //----------------------------------------------------------------------------------------------
515 Double_t PndDrcHitProducerReal::FindPhiRot(Double_t xx, Double_t yy){ // returns [degrees]
516 
517  TVector3 hit;
518  hit.SetXYZ(xx,yy,0.);
519  Double_t startPhi = hit.Phi()/fpi*180.; // [degrees]
520  if(startPhi < 0.){startPhi = 360. + hit.Phi()*180./fpi;}
521  //cout<<"-I- FindPhoRot: start phi = "<<startPhi<<endl;
522  //cout<<"-I- InBarCoordinateSystem: dphi = "<<fDphi<<endl;
523  Double_t PhiRot = 0.; //[degrees]
524  if(startPhi >= 0. && startPhi < 90.){
525  PhiRot = TMath::Floor(startPhi/fdphi) *fdphi + fdphi/2.;
526  }
527  if(startPhi >= 90. && startPhi < 270.){
528  PhiRot = 90. + fpipehAngle + TMath::Floor((startPhi-90.-fpipehAngle)/fdphi) *fdphi + fdphi/2.;
529  }
530  if(startPhi >= 270. && startPhi < 360.){
531  PhiRot = 270. + fpipehAngle + TMath::Floor((startPhi-270.-fpipehAngle)/fdphi) *fdphi + fdphi/2.;
532  }
533  //cout<<"-I- FindPhiRot: PhiRot = "<<PhiRot<<endl;
534  return PhiRot; // degrees
535 }
536 
537 //----------------------------------------------------------------------------------------------------------
538 Double_t PndDrcHitProducerReal::FindOutPoint(Double_t x0, Double_t xEn, Double_t a, Double_t *NN, Bool_t print){
539  Double_t m=99.;
540  Double_t n=TMath::Floor(x0/a);
541  m = n;
542  if(print){std::cout<<"n = "<<n<<", NN = "<<*NN<<", x0 = "<<x0<<", a = "<<a<<std::endl;}
543  Double_t x1 = x0 - n*a;
544  if(x0 < 0.){x1 = x0 - (n+1)*a;}
545  if(print){std::cout<<"xy = "<< x1<<std::endl;}
546  Double_t xK = 0.;
547  if((m/2. - TMath::Floor(m/2.)) == 0.) { // 4etnoe
548  if(print){std::cout<<"odd==0"<<std::endl;}
549  if(x0 >= 0. && x1 + xEn <= a){xK = x1 + xEn;}
550  if(x0 >= 0. && x1 + xEn > a){xK = 2*a - x1 - xEn; n = 1. + n;}
551  if(x0 < 0. && x1 + xEn >= 0.){xK = a - (x1 + xEn); n = -1. -n;}
552  if(x0 < 0. && x1 + xEn < 0.) {xK = a + x1 + xEn; n = -n;}
553  if(print){std::cout<<"xK = "<< xK<<", n = "<<n<<std::endl;}
554 
555  }
556 
557  if((m/2. - TMath::Floor(m/2.)) != 0.) { // ne4etnoe
558  if(print){std::cout<<"even!=0"<<std::endl;}
559  if(x0 >= 0. && x1 + xEn <= a){xK = a - (x1 + xEn);}
560  if(x0 >= 0. && x1 + xEn > a){xK = x1 + xEn - a; n = 1. + n;}
561  if(x0 < 0. && x1 + xEn >= 0.){xK = x1 + xEn; n = -1. -n;}
562  if(x0 < 0. && x1 + xEn < 0.) {xK = - (x1 + xEn); n = -n;}
563  if(print){std::cout<<"xK = "<< xK<<", n = "<<n<<std::endl;}
564 
565  }
566 
567  *NN = n;
568  return xK;
569 }
570 
571 // ----- Add Hit to HitCollection --------------------------------------
572 PndDrcHit* PndDrcHitProducerReal::AddHit(Int_t detID,
573  TVector3 posHit,
574  TVector3 dPosHit,
575  Double_t thetaC,
576  Double_t errThetaC,
577  Int_t index){
578  TClonesArray& clref = *fHitArray;
579  Int_t size = clref.GetEntriesFast();
580  return new(clref[size]) PndDrcHit(detID, detID,
581  posHit,
582  dPosHit,
583  thetaC,
584  errThetaC,
585  index);
586 }
587 // ----- Add Photon Detector Hit to HitCollection --------------------------------------
588 PndDrcPDHit* PndDrcHitProducerReal::AddPDHit(Int_t detID,
589  TVector3 posHit,
590  TVector3 dPosHit,
591  Double_t time,
592  Double_t timeThreshold,
593  Int_t index1){
594  TClonesArray& clrefPD = *fPDHitArray;
595  Int_t size = clrefPD.GetEntriesFast();
596  return new(clrefPD[size]) PndDrcPDHit(detID,
597  detID,
598  posHit,
599  dPosHit,
600  time,
601  timeThreshold,
602  index1);
603 }
604 // ----- Set Photon Detector parameter -------------------------------------------
605 void PndDrcHitProducerReal::SetPhotonDetEff(Double_t& fLambdaMin,
606  Double_t& fLambdaMax, Double_t& fLambdaStep, Double_t fEfficiency[])
607 {
608  // if (fVerbose > 0) cout << "SetPhotoDetParamter called for Photocathode type " << fDetType << endl;
609 
610  if (fDetType == 1){
611 
614  fLambdaMin = 300.;
615  fLambdaMax = 660.;
616  fLambdaStep = 10.;
617 
618  fEfficiency[0] = 0.1;
619  fEfficiency[1] = 0.12;
620  fEfficiency[2] = 0.13;
621  fEfficiency[3] = 0.15;
622  fEfficiency[4] = 0.18;
623  fEfficiency[5] = 0.2;
624  fEfficiency[6] = 0.22;
625  fEfficiency[7] = 0.24;
626  fEfficiency[8] = 0.26;
627  fEfficiency[9] = 0.28;
628  fEfficiency[10] = 0.29;
629  fEfficiency[11] = 0.288;
630  fEfficiency[12] = 0.286;
631  fEfficiency[13] = 0.284;
632  fEfficiency[14] = 0.282;
633  fEfficiency[15] = 0.28;
634  fEfficiency[16] = 0.26;
635  fEfficiency[17] = 0.24;
636  fEfficiency[18] = 0.2;
637  fEfficiency[19] = 0.17;
638  fEfficiency[20] = 0.15;
639  fEfficiency[21] = 0.13;
640  fEfficiency[22] = 0.1;
641  fEfficiency[23] = 0.09;
642  fEfficiency[24] = 0.07;
643  fEfficiency[25] = 0.06;
644  fEfficiency[26] = 0.05;
645  fEfficiency[27] = 0.04;
646  fEfficiency[28] = 0.03;
647  fEfficiency[29] = 0.02;
648  fEfficiency[30] = 0.015;
649  fEfficiency[31] = 0.01;
650  fEfficiency[32] = 0.005;
651  fEfficiency[33] = 0.002;
652  fEfficiency[34] = 0.001;
653  fEfficiency[35] = 0.;
654  }
655  else {
656  cout << "ERROR: photocathode type not specified" << endl;
657 
658  fLambdaMin = 100.;
659  fLambdaMax = 100.;
660  fLambdaStep = 100.;
661 
662  fEfficiency[0] = 0.;
663 
664  }
665 
666 }
667 //-------------------------------------------------------------------------------
668 void PndDrcHitProducerReal::SetFakeDetEff(Double_t& fLambdaMin,
669  Double_t& fLambdaMax, Double_t& fLambdaStep, Double_t fEfficiency[])
670 {
671  // if (fVerbose > 0) cout << "SetPhotoDetParamter called for Photocathode type " << fDetType << endl;
672 
673  if (fDetType == 1){
674  fLambdaMin = 300.;
675  fLambdaMax = 700.;
676  fLambdaStep = 400.;
677 
678  fEfficiency[0] = 1.;
679  fEfficiency[1] = 1.;
680  }else{
681  cout << "ERROR: photocathode type not specified" << endl;
682 
683  fLambdaMin = 300.;
684  fLambdaMax = 700.;
685  fLambdaStep = 400.;
686 
687  fEfficiency[0] = 0.;
688  }
689 }
690 
691 void PndDrcHitProducerReal::SetPhotonDetEffNew(Double_t& fLambdaMin,
692  Double_t& fLambdaMax, Double_t& fLambdaStep, Double_t fEfficiency[])
693 {
694  // if (fVerbose > 0) cout << "SetPhotoDetParamter called for Photocathode type " << fDetType << endl;
695 
696  if (fDetType == 1){
697 
698  fLambdaMin = 200.;
699  fLambdaMax = 700.;
700  fLambdaStep = 1.;
701 
702  Float_t credibleLimit=280.;
703 
704  // quantum efficiency data from Alex Britting, Jan 25, 2011
705  // unit is percent
706  // first value is at 200 nm, last at 700 nm
707  // credible range start around 250nm, >= 280nm to be safe
708 
709  fEfficiency[0]= 231.84;
710  fEfficiency[1]= 615.36;
711  fEfficiency[2]= 657.4;
712  fEfficiency[3]= 258.78;
713  fEfficiency[4]= 9839.92;
714  fEfficiency[5]= 44.67;
715  fEfficiency[6]= 67.87;
716  fEfficiency[7]= 51.01;
717  fEfficiency[8]= 41.49;
718  fEfficiency[9]= 5.36;
719  fEfficiency[10]= 49.4;
720  fEfficiency[11]= 2.13;
721  fEfficiency[12]= 35.49;
722  fEfficiency[13]= 8.66;
723  fEfficiency[14]= 5.03;
724  fEfficiency[15]= 7.51;
725  fEfficiency[16]= 13.27;
726  fEfficiency[17]= 18.71;
727  fEfficiency[18]= 3.92;
728  fEfficiency[19]= 3.66;
729  fEfficiency[20]= 8.2;
730  fEfficiency[21]= 0.56;
731  fEfficiency[22]= 7.68;
732  fEfficiency[23]= 2.87;
733  fEfficiency[24]= 10.06;
734  fEfficiency[25]= 3.47;
735  fEfficiency[26]= 3.39;
736  fEfficiency[27]= 6.99;
737  fEfficiency[28]= 6.01;
738  fEfficiency[29]= 4.92;
739  fEfficiency[30]= 6.25;
740  fEfficiency[31]= 5.97;
741  fEfficiency[32]= 6.92;
742  fEfficiency[33]= 8.29;
743  fEfficiency[34]= 10.45;
744  fEfficiency[35]= 8.68;
745  fEfficiency[36]= 8.6;
746  fEfficiency[37]= 9.79;
747  fEfficiency[38]= 11.76;
748  fEfficiency[39]= 9.53;
749  fEfficiency[40]= 10.98;
750  fEfficiency[41]= 9.9;
751  fEfficiency[42]= 10.97;
752  fEfficiency[43]= 11.31;
753  fEfficiency[44]= 10.88;
754  fEfficiency[45]= 10.78;
755  fEfficiency[46]= 12.16;
756  fEfficiency[47]= 12.38;
757  fEfficiency[48]= 12.37;
758  fEfficiency[49]= 13.04;
759  fEfficiency[50]= 12.36;
760  fEfficiency[51]= 13.18;
761  fEfficiency[52]= 13.7;
762  fEfficiency[53]= 13.85;
763  fEfficiency[54]= 13.66;
764  fEfficiency[55]= 13.98;
765  fEfficiency[56]= 14.55;
766  fEfficiency[57]= 14.93;
767  fEfficiency[58]= 14.82;
768  fEfficiency[59]= 14.97;
769  fEfficiency[60]= 14.98;
770  fEfficiency[61]= 15.14;
771  fEfficiency[62]= 15.35;
772  fEfficiency[63]= 15.37;
773  fEfficiency[64]= 15.43;
774  fEfficiency[65]= 15.49;
775  fEfficiency[66]= 15.59;
776  fEfficiency[67]= 15.84;
777  fEfficiency[68]= 15.84;
778  fEfficiency[69]= 15.92;
779  fEfficiency[70]= 16.01;
780  fEfficiency[71]= 16.22;
781  fEfficiency[72]= 16.41;
782  fEfficiency[73]= 16.42;
783  fEfficiency[74]= 16.52;
784  fEfficiency[75]= 16.86;
785  fEfficiency[76]= 17.1;
786  fEfficiency[77]= 17.17;
787  fEfficiency[78]= 17.22;
788  fEfficiency[79]= 17.46;
789  fEfficiency[80]= 17.79;
790  fEfficiency[81]= 17.99;
791  fEfficiency[82]= 18.13;
792  fEfficiency[83]= 18.33;
793  fEfficiency[84]= 18.34;
794  fEfficiency[85]= 18.53;
795  fEfficiency[86]= 18.72;
796  fEfficiency[87]= 18.95;
797  fEfficiency[88]= 19.02;
798  fEfficiency[89]= 19.15;
799  fEfficiency[90]= 19.28;
800  fEfficiency[91]= 19.45;
801  fEfficiency[92]= 19.66;
802  fEfficiency[93]= 19.69;
803  fEfficiency[94]= 19.77;
804  fEfficiency[95]= 19.73;
805  fEfficiency[96]= 19.95;
806  fEfficiency[97]= 19.98;
807  fEfficiency[98]= 20.17;
808  fEfficiency[99]= 20.29;
809  fEfficiency[100]=20.33;
810  fEfficiency[101]=20.37;
811  fEfficiency[102]=20.47;
812  fEfficiency[103]=20.48;
813  fEfficiency[104]=20.57;
814  fEfficiency[105]=20.75;
815  fEfficiency[106]=20.8;
816  fEfficiency[107]=20.84;
817  fEfficiency[108]=20.86;
818  fEfficiency[109]=20.88;
819  fEfficiency[110]=21.0;
820  fEfficiency[111]=21.06;
821  fEfficiency[112]=21.0;
822  fEfficiency[113]=21.06;
823  fEfficiency[114]=21.06;
824  fEfficiency[115]=21.04;
825  fEfficiency[116]=21.1;
826  fEfficiency[117]=21.14;
827  fEfficiency[118]=21.08;
828  fEfficiency[119]=21.17;
829  fEfficiency[120]=21.3;
830  fEfficiency[121]=21.38;
831  fEfficiency[122]=21.49;
832  fEfficiency[123]=21.58;
833  fEfficiency[124]=21.69;
834  fEfficiency[125]=21.77;
835  fEfficiency[126]=21.87;
836  fEfficiency[127]=22.02;
837  fEfficiency[128]=22.13;
838  fEfficiency[129]=22.29;
839  fEfficiency[130]=22.35;
840  fEfficiency[131]=22.45;
841  fEfficiency[132]=22.53;
842  fEfficiency[133]=22.55;
843  fEfficiency[134]=22.64;
844  fEfficiency[135]=22.67;
845  fEfficiency[136]=22.73;
846  fEfficiency[137]=22.74;
847  fEfficiency[138]=22.71;
848  fEfficiency[139]=22.79;
849  fEfficiency[140]=22.76;
850  fEfficiency[141]=22.77;
851  fEfficiency[142]=22.76;
852  fEfficiency[143]=22.75;
853  fEfficiency[144]=22.78;
854  fEfficiency[145]=22.7;
855  fEfficiency[146]=22.68;
856  fEfficiency[147]=22.72;
857  fEfficiency[148]=22.66;
858  fEfficiency[149]=22.64;
859  fEfficiency[150]=22.7;
860  fEfficiency[151]=22.67;
861  fEfficiency[152]=22.71;
862  fEfficiency[153]=22.67;
863  fEfficiency[154]=22.75;
864  fEfficiency[155]=22.77;
865  fEfficiency[156]=22.83;
866  fEfficiency[157]=22.84;
867  fEfficiency[158]=22.93;
868  fEfficiency[159]=22.97;
869  fEfficiency[160]=23.0;
870  fEfficiency[161]=23.08;
871  fEfficiency[162]=23.16;
872  fEfficiency[163]=23.27;
873  fEfficiency[164]=23.25;
874  fEfficiency[165]=23.37;
875  fEfficiency[166]=23.44;
876  fEfficiency[167]=23.49;
877  fEfficiency[168]=23.55;
878  fEfficiency[169]=23.52;
879  fEfficiency[170]=23.58;
880  fEfficiency[171]=23.64;
881  fEfficiency[172]=23.63;
882  fEfficiency[173]=23.58;
883  fEfficiency[174]=23.64;
884  fEfficiency[175]=23.63;
885  fEfficiency[176]=23.62;
886  fEfficiency[177]=23.64;
887  fEfficiency[178]=23.63;
888  fEfficiency[179]=23.66;
889  fEfficiency[180]=23.59;
890  fEfficiency[181]=23.59;
891  fEfficiency[182]=23.56;
892  fEfficiency[183]=23.58;
893  fEfficiency[184]=23.63;
894  fEfficiency[185]=23.57;
895  fEfficiency[186]=23.66;
896  fEfficiency[187]=23.62;
897  fEfficiency[188]=23.67;
898  fEfficiency[189]=23.64;
899  fEfficiency[190]=23.54;
900  fEfficiency[191]=23.57;
901  fEfficiency[192]=23.51;
902  fEfficiency[193]=23.53;
903  fEfficiency[194]=23.45;
904  fEfficiency[195]=23.3;
905  fEfficiency[196]=23.41;
906  fEfficiency[197]=23.25;
907  fEfficiency[198]=23.21;
908  fEfficiency[199]=23.08;
909  fEfficiency[200]=23.01;
910  fEfficiency[201]=22.92;
911  fEfficiency[202]=22.9;
912  fEfficiency[203]=22.76;
913  fEfficiency[204]=22.76;
914  fEfficiency[205]=22.61;
915  fEfficiency[206]=22.53;
916  fEfficiency[207]=22.48;
917  fEfficiency[208]=22.39;
918  fEfficiency[209]=22.29;
919  fEfficiency[210]=22.24;
920  fEfficiency[211]=22.2;
921  fEfficiency[212]=22.12;
922  fEfficiency[213]=22.07;
923  fEfficiency[214]=21.96;
924  fEfficiency[215]=21.89;
925  fEfficiency[216]=21.87;
926  fEfficiency[217]=21.76;
927  fEfficiency[218]=21.74;
928  fEfficiency[219]=21.58;
929  fEfficiency[220]=21.49;
930  fEfficiency[221]=21.48;
931  fEfficiency[222]=21.37;
932  fEfficiency[223]=21.29;
933  fEfficiency[224]=21.2;
934  fEfficiency[225]=21.17;
935  fEfficiency[226]=21.03;
936  fEfficiency[227]=20.98;
937  fEfficiency[228]=20.92;
938  fEfficiency[229]=20.85;
939  fEfficiency[230]=20.76;
940  fEfficiency[231]=20.69;
941  fEfficiency[232]=20.58;
942  fEfficiency[233]=20.56;
943  fEfficiency[234]=20.47;
944  fEfficiency[235]=20.37;
945  fEfficiency[236]=20.32;
946  fEfficiency[237]=20.24;
947  fEfficiency[238]=20.13;
948  fEfficiency[239]=20.08;
949  fEfficiency[240]=19.9;
950  fEfficiency[241]=19.84;
951  fEfficiency[242]=19.77;
952  fEfficiency[243]=19.69;
953  fEfficiency[244]=19.63;
954  fEfficiency[245]=19.51;
955  fEfficiency[246]=19.41;
956  fEfficiency[247]=19.27;
957  fEfficiency[248]=19.06;
958  fEfficiency[249]=19.01;
959  fEfficiency[250]=18.87;
960  fEfficiency[251]=18.7;
961  fEfficiency[252]=18.49;
962  fEfficiency[253]=18.41;
963  fEfficiency[254]=18.17;
964  fEfficiency[255]=17.98;
965  fEfficiency[256]=17.84;
966  fEfficiency[257]=17.69;
967  fEfficiency[258]=17.5;
968  fEfficiency[259]=17.25;
969  fEfficiency[260]=17.15;
970  fEfficiency[261]=16.98;
971  fEfficiency[262]=16.79;
972  fEfficiency[263]=16.66;
973  fEfficiency[264]=16.48;
974  fEfficiency[265]=16.32;
975  fEfficiency[266]=16.19;
976  fEfficiency[267]=16.02;
977  fEfficiency[268]=15.88;
978  fEfficiency[269]=15.77;
979  fEfficiency[270]=15.67;
980  fEfficiency[271]=15.5;
981  fEfficiency[272]=15.39;
982  fEfficiency[273]=15.23;
983  fEfficiency[274]=15.09;
984  fEfficiency[275]=15.04;
985  fEfficiency[276]=14.92;
986  fEfficiency[277]=14.75;
987  fEfficiency[278]=14.7;
988  fEfficiency[279]=14.5;
989  fEfficiency[280]=14.45;
990  fEfficiency[281]=14.34;
991  fEfficiency[282]=14.25;
992  fEfficiency[283]=14.16;
993  fEfficiency[284]=14.13;
994  fEfficiency[285]=14.0;
995  fEfficiency[286]=13.92;
996  fEfficiency[287]=13.84;
997  fEfficiency[288]=13.76;
998  fEfficiency[289]=13.73;
999  fEfficiency[290]=13.61;
1000  fEfficiency[291]=13.54;
1001  fEfficiency[292]=13.52;
1002  fEfficiency[293]=13.45;
1003  fEfficiency[294]=13.41;
1004  fEfficiency[295]=13.39;
1005  fEfficiency[296]=13.31;
1006  fEfficiency[297]=13.22;
1007  fEfficiency[298]=13.17;
1008  fEfficiency[299]=13.13;
1009  fEfficiency[300]=13.06;
1010  fEfficiency[301]=13.2;
1011  fEfficiency[302]=13.09;
1012  fEfficiency[303]=12.97;
1013  fEfficiency[304]=12.92;
1014  fEfficiency[305]=12.73;
1015  fEfficiency[306]=12.65;
1016  fEfficiency[307]=12.4;
1017  fEfficiency[308]=12.22;
1018  fEfficiency[309]=12.02;
1019  fEfficiency[310]=11.79;
1020  fEfficiency[311]=11.59;
1021  fEfficiency[312]=11.33;
1022  fEfficiency[313]=11.03;
1023  fEfficiency[314]=10.68;
1024  fEfficiency[315]=10.46;
1025  fEfficiency[316]=10.14;
1026  fEfficiency[317]=9.88;
1027  fEfficiency[318]=9.62;
1028  fEfficiency[319]=9.36;
1029  fEfficiency[320]=9.14;
1030  fEfficiency[321]=8.87;
1031  fEfficiency[322]=8.63;
1032  fEfficiency[323]=8.51;
1033  fEfficiency[324]=8.24;
1034  fEfficiency[325]=8.07;
1035  fEfficiency[326]=7.88;
1036  fEfficiency[327]=7.77;
1037  fEfficiency[328]=7.65;
1038  fEfficiency[329]=7.52;
1039  fEfficiency[330]=7.35;
1040  fEfficiency[331]=7.27;
1041  fEfficiency[332]=7.21;
1042  fEfficiency[333]=7.1;
1043  fEfficiency[334]=6.92;
1044  fEfficiency[335]=6.89;
1045  fEfficiency[336]=6.79;
1046  fEfficiency[337]=6.74;
1047  fEfficiency[338]=6.56;
1048  fEfficiency[339]=6.54;
1049  fEfficiency[340]=6.5;
1050  fEfficiency[341]=6.39;
1051  fEfficiency[342]=6.33;
1052  fEfficiency[343]=6.25;
1053  fEfficiency[344]=6.27;
1054  fEfficiency[345]=6.14;
1055  fEfficiency[346]=6.06;
1056  fEfficiency[347]=6.04;
1057  fEfficiency[348]=6.01;
1058  fEfficiency[349]=5.91;
1059  fEfficiency[350]=5.89;
1060  fEfficiency[351]=5.79;
1061  fEfficiency[352]=5.75;
1062  fEfficiency[353]=5.75;
1063  fEfficiency[354]=5.67;
1064  fEfficiency[355]=5.61;
1065  fEfficiency[356]=5.51;
1066  fEfficiency[357]=5.52;
1067  fEfficiency[358]=5.43;
1068  fEfficiency[359]=5.43;
1069  fEfficiency[360]=5.34;
1070  fEfficiency[361]=5.31;
1071  fEfficiency[362]=5.35;
1072  fEfficiency[363]=5.23;
1073  fEfficiency[364]=5.2;
1074  fEfficiency[365]=5.14;
1075  fEfficiency[366]=5.11;
1076  fEfficiency[367]=5.11;
1077  fEfficiency[368]=5.01;
1078  fEfficiency[369]=4.98;
1079  fEfficiency[370]=4.93;
1080  fEfficiency[371]=4.99;
1081  fEfficiency[372]=4.89;
1082  fEfficiency[373]=4.82;
1083  fEfficiency[374]=4.87;
1084  fEfficiency[375]=4.8;
1085  fEfficiency[376]=4.7;
1086  fEfficiency[377]=4.65;
1087  fEfficiency[378]=4.65;
1088  fEfficiency[379]=4.61;
1089  fEfficiency[380]=4.49;
1090  fEfficiency[381]=4.56;
1091  fEfficiency[382]=4.44;
1092  fEfficiency[383]=4.42;
1093  fEfficiency[384]=4.44;
1094  fEfficiency[385]=4.35;
1095  fEfficiency[386]=4.35;
1096  fEfficiency[387]=4.27;
1097  fEfficiency[388]=4.29;
1098  fEfficiency[389]=4.19;
1099  fEfficiency[390]=4.13;
1100  fEfficiency[391]=4.08;
1101  fEfficiency[392]=4.02;
1102  fEfficiency[393]=4.07;
1103  fEfficiency[394]=3.92;
1104  fEfficiency[395]=3.95;
1105  fEfficiency[396]=3.88;
1106  fEfficiency[397]=3.82;
1107  fEfficiency[398]=3.86;
1108  fEfficiency[399]=3.74;
1109  fEfficiency[400]=3.71;
1110  fEfficiency[401]=3.66;
1111  fEfficiency[402]=3.72;
1112  fEfficiency[403]=3.62;
1113  fEfficiency[404]=3.55;
1114  fEfficiency[405]=3.56;
1115  fEfficiency[406]=3.57;
1116  fEfficiency[407]=3.45;
1117  fEfficiency[408]=3.38;
1118  fEfficiency[409]=3.36;
1119  fEfficiency[410]=3.36;
1120  fEfficiency[411]=3.28;
1121  fEfficiency[412]=3.25;
1122  fEfficiency[413]=3.19;
1123  fEfficiency[414]=3.26;
1124  fEfficiency[415]=3.13;
1125  fEfficiency[416]=3.17;
1126  fEfficiency[417]=3.15;
1127  fEfficiency[418]=3.04;
1128  fEfficiency[419]=2.98;
1129  fEfficiency[420]=2.93;
1130  fEfficiency[421]=2.98;
1131  fEfficiency[422]=2.9;
1132  fEfficiency[423]=2.89;
1133  fEfficiency[424]=2.9;
1134  fEfficiency[425]=2.81;
1135  fEfficiency[426]=2.74;
1136  fEfficiency[427]=2.81;
1137  fEfficiency[428]=2.68;
1138  fEfficiency[429]=2.73;
1139  fEfficiency[430]=2.7;
1140  fEfficiency[431]=2.57;
1141  fEfficiency[432]=2.58;
1142  fEfficiency[433]=2.55;
1143  fEfficiency[434]=2.55;
1144  fEfficiency[435]=2.37;
1145  fEfficiency[436]=2.39;
1146  fEfficiency[437]=2.39;
1147  fEfficiency[438]=2.44;
1148  fEfficiency[439]=2.37;
1149  fEfficiency[440]=2.26;
1150  fEfficiency[441]=2.27;
1151  fEfficiency[442]=2.27;
1152  fEfficiency[443]=2.23;
1153  fEfficiency[444]=2.26;
1154  fEfficiency[445]=2.14;
1155  fEfficiency[446]=2.08;
1156  fEfficiency[447]=2.15;
1157  fEfficiency[448]=2.06;
1158  fEfficiency[449]=2.09;
1159  fEfficiency[450]=2.04;
1160  fEfficiency[451]=2.0;
1161  fEfficiency[452]=1.95;
1162  fEfficiency[453]=2.02;
1163  fEfficiency[454]=1.87;
1164  fEfficiency[455]=1.9;
1165  fEfficiency[456]=1.8;
1166  fEfficiency[457]=1.87;
1167  fEfficiency[458]=1.85;
1168  fEfficiency[459]=1.87;
1169  fEfficiency[460]=1.81;
1170  fEfficiency[461]=1.86;
1171  fEfficiency[462]=1.74;
1172  fEfficiency[463]=1.74;
1173  fEfficiency[464]=1.63;
1174  fEfficiency[465]=1.59;
1175  fEfficiency[466]=1.5;
1176  fEfficiency[467]=1.5;
1177  fEfficiency[468]=1.44;
1178  fEfficiency[469]=1.47;
1179  fEfficiency[470]=1.32;
1180  fEfficiency[471]=1.24;
1181  fEfficiency[472]=1.28;
1182  fEfficiency[473]=1.19;
1183  fEfficiency[474]=1.21;
1184  fEfficiency[475]=1.21;
1185  fEfficiency[476]=1.1;
1186  fEfficiency[477]=1.1;
1187  fEfficiency[478]=1.05;
1188  fEfficiency[479]=1.06;
1189  fEfficiency[480]=0.94;
1190  fEfficiency[481]=0.92;
1191  fEfficiency[482]=0.87;
1192  fEfficiency[483]=0.92;
1193  fEfficiency[484]=0.81;
1194  fEfficiency[485]=0.86;
1195  fEfficiency[486]=0.78;
1196  fEfficiency[487]=0.77;
1197  fEfficiency[488]=0.8;
1198  fEfficiency[489]=0.67;
1199  fEfficiency[490]=0.7;
1200  fEfficiency[491]=0.81;
1201  fEfficiency[492]=0.61;
1202  fEfficiency[493]=0.64;
1203  fEfficiency[494]=0.71;
1204  fEfficiency[495]=0.66;
1205  fEfficiency[496]=0.67;
1206  fEfficiency[497]=0.68;
1207  fEfficiency[498]=0.69;
1208  fEfficiency[499]=0.68;
1209  fEfficiency[500]=0.73;
1210 
1211  // still need to convert from percent and cut values below credible limit
1212  for (Int_t iBin=0;iBin<501;iBin++)
1213  {
1214  if (iBin<(Int_t)(credibleLimit-fLambdaMin))
1215  {
1216  fEfficiency[iBin]=0.;
1217  }
1218  else
1219  {
1220  fEfficiency[iBin]=fEfficiency[iBin]/100.;
1221  }
1222  // cout << "efficiency is " << fEfficiency[iBin] << " at " << fLambdaMin+iBin<<endl;
1223  }
1224 
1225  }
1226  else {
1227  cout << "ERROR: photocathode type not specified" << endl;
1228 
1229  fLambdaMin = 100.;
1230  fLambdaMax = 100.;
1231  fLambdaStep = 100.;
1232 
1233  fEfficiency[0] = 0.;
1234 
1235  }
1236 
1237 }
1238 //-------------------------------------------------------------------------------------------
1239 void PndDrcHitProducerReal::SetPhotonDetEffOld(Double_t& fLambdaMin,
1240  Double_t& fLambdaMax, Double_t& fLambdaStep, Double_t fEfficiency[])
1241 {
1242  // if (fVerbose > 0) cout << "SetPhotoDetParamter called for Photocathode type " << fDetType << endl;
1243 
1244  if (fDetType == 1){
1245 
1249  fLambdaMin = 280.;
1250  fLambdaMax = 600.;
1251  fLambdaStep = 10.;
1252 
1253  fEfficiency[0] = 0.020;
1254  fEfficiency[1] = 0.050;
1255  fEfficiency[2] = 0.110;
1256  fEfficiency[3] = 0.170;
1257  fEfficiency[4] = 0.225;
1258  fEfficiency[5] = 0.260;
1259  fEfficiency[6] = 0.275;
1260  fEfficiency[7] = 0.290;
1261  fEfficiency[8] = 0.310;
1262  fEfficiency[9] = 0.305;
1263  fEfficiency[10] = 0.302;
1264  fEfficiency[11] = 0.290;
1265  fEfficiency[12] = 0.279;
1266  fEfficiency[13] = 0.268;
1267  fEfficiency[14] = 0.262;
1268  fEfficiency[15] = 0.250;
1269  fEfficiency[16] = 0.240;
1270  fEfficiency[17] = 0.225;
1271  fEfficiency[18] = 0.210;
1272  fEfficiency[19] = 0.195;
1273  fEfficiency[20] = 0.178;
1274  fEfficiency[21] = 0.170;
1275  fEfficiency[22] = 0.155;
1276  fEfficiency[23] = 0.130;
1277  fEfficiency[24] = 0.112;
1278  fEfficiency[25] = 0.085;
1279  fEfficiency[26] = 0.066;
1280  fEfficiency[27] = 0.055;
1281  fEfficiency[28] = 0.040;
1282  fEfficiency[29] = 0.030;
1283  fEfficiency[30] = 0.025;
1284  fEfficiency[31] = 0.018;
1285  fEfficiency[32] = 0.010;
1286  }
1287  else {
1288  cout << "ERROR: photocathode type not specified" << endl;
1289 
1290  fLambdaMin = 100.;
1291  fLambdaMax = 100.;
1292  fLambdaStep = 100.;
1293 
1294  fEfficiency[0] = 0.;
1295 
1296  }
1297 
1298 }
1299 //------------------------------------------------------------------------------------
1300 //void PndDrcHitProducerReal::SetPhotonTransportEff(Double_t lambda, Double_t angle){
1301 void PndDrcHitProducerReal::SetPhotonTransportEff(Double_t& fLambdaMin,
1302  Double_t& fLambdaMax, Double_t& fLambdaStep, Double_t& fAngleStep, Int_t& fLambdaPoints, Double_t fEfficiency[])
1303 {
1304  if (fDetType == 1){
1305 
1306  //cout<<"-I- SetPhotonTransportEff: lam = "<<lambda<<", angle = "<<angle<<endl;
1307 
1308  fLambdaMin = 280.;
1309  fLambdaMax = 650.;
1310  fLambdaStep = 10.;
1311 
1312  fAngleStep = fpi/2./20.;
1313  fLambdaPoints = 38;
1314 
1315  // refraction index of silica:
1316  TF1 *d1 = new TF1("d1", "sqrt(1 + ([0]*x^2/(x^2-[1]^2)) + ([2]*x^2/(x^2-[3]^2)) + ([4]*x^2/(x^2-[5]^2)))",fLambdaMin/1000.,fLambdaMax/1000.);
1317  d1->SetParameters(0.696, 0.068, 0.407, 0.116, 0.897, 9.896);
1318 
1319  // reflection probability according to the scalar theory
1320  TF2* d3 = new TF2("d3", "1. - pow(4.*3.14159*cos(y)*[0]*d1/x,2)",fLambdaMin/1000.,fLambdaMax/1000., 0.,fpi/2.);
1321  d3->SetParameter(0,fRoughness);
1322 
1323  //cout<<"-I- SetPhotonTransportEff: reflection coef(l = 0.53, ang = 0.95) = "<<d3->Eval(0.53,0.95)<<endl;
1324 
1325  for(int iang=0; iang<21; iang++){ // lambda range, 38 points, step = (0.65-0.28)/37 = 0.01 = 10 nm
1326  for(int ilam=0; ilam< 38; ilam++){ // angle range, 21 points, step = pi/2/20
1327  fEfficiency[iang * 38 + ilam] = d3->Eval(fLambdaMin/1000.+ilam*0.01, iang*fpi/2./20.);
1328  }
1329  }
1330 
1331  delete(d1);
1332  delete(d3);
1333  //cout<<"-I- SetPhotonTransportEff: reflection coef = "<<d3->Eval(lambda, angle)<<endl;
1334 
1335  //d3->Eval(lambda, angle);
1336 
1337  }
1338  else {
1339  cout << "ERROR: photocathode type not specified" << endl;
1340 
1341  }
1342 
1343 }
1344 // ----- Find Photon Hit Position---------------------------------------------------
1345 void PndDrcHitProducerReal::FindDrcHitPosition(Double_t xPoint, Double_t yPoint,
1346  Double_t& xHit, Double_t& yHit, Int_t& ) // pmtID //[R.K.03/2017] unused variable(s)
1347 {
1348  Double_t pixelDim = fPixelDim;
1349  if(xPoint >= 0)
1350  xHit= pixelDim/2. + pixelDim*((Int_t)(xPoint/pixelDim));
1351  else
1352  xHit= -pixelDim/2. + pixelDim*((Int_t)(xPoint/pixelDim));
1353 
1354  if(yPoint >= 0)
1355  yHit= pixelDim/2. + pixelDim*((Int_t)(yPoint/pixelDim));
1356  else
1357  yHit= -pixelDim/2. + pixelDim*((Int_t)(yPoint/pixelDim));
1358 
1359 }
1360 
1361 // ----- Find Photon Hit Position Tilt---------------------------------------------------
1362 void PndDrcHitProducerReal::FindDrcHitPositionTilt(Double_t xPoint, Double_t yPoint,
1363  Double_t& xHit, Double_t& yHit, Int_t& ) // pmtID //[R.K.03/2017] unused variable(s)
1364 {
1365  Double_t pixelDim = fPixelDim;
1366  Double_t pixelDimY = fPixelDim*cos(ftilt/180.*fpi);
1367  //cout<<"-I- HitProducerReal: ftilt = "<<ftilt<<", pixelDimY = "<<pixelDimY<<endl;
1368  if(xPoint >= 0)
1369  xHit= pixelDim/2. + pixelDim*((Int_t)(xPoint/pixelDim));
1370  else
1371  xHit= -pixelDim/2. + pixelDim*((Int_t)(xPoint/pixelDim));
1372 
1373  if(yPoint >= 0)
1374  yHit= pixelDimY/2. + pixelDimY*((Int_t)(yPoint/pixelDimY));
1375  else
1376  yHit= -pixelDimY/2. + pixelDimY*((Int_t)(yPoint/pixelDimY));
1377 
1378 }
1379 
1380 //-------------Smear Time------------------------------------
1381 void PndDrcHitProducerReal::Smear(Double_t& time, Double_t sigt){
1382  Double_t t=time;
1383  Double_t dt;
1384  dt=gRandom->Gaus(0,sigt);
1385  t +=dt;
1386  time=t;
1387  return;
1388 }
1389 
1390 // ----- Finish Task ---------------------------------------------------
1391 void PndDrcHitProducerReal::Finish()
1392 {
1393  cout << "-I- PndDrcHitProducerReal: Finish" << endl;
1394 }
1395 // -------------------------------------------------------------------------
1396 
1397 
1398 ClassImp(PndDrcHitProducerReal)
x0
Double_t x0
Definition: checkhelixhit.C:70
PndDrcBarPoint::GetMass
Double_t GetMass() const
Definition: PndDrcBarPoint.h:45
PndGeoDrc::nEV
Double_t nEV()
Definition: PndGeoDrc.h:76
fVerbose
int fVerbose
Definition: poormantracks.C:24
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
rand
TRandom1 rand
Definition: RiemannTest.C:20
PndDrcHitProducerReal::SetPhotonDetEff
void SetPhotonDetEff(Double_t &lambda_min, Double_t &lambda_max, Double_t &lambda_step, Double_t efficiency[])
Definition: PndDrcHitProducerReal.cxx:605
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: P4_F32vec4.h:109
PndGeoDrc::BBoxNum
Double_t BBoxNum()
Definition: PndGeoDrc.h:136
m
__m128 m
Definition: P4_F32vec4.h:28
PndDrcHitProducerReal::Smear
void Smear(Double_t &time, Double_t sigt)
Definition: PndDrcHitProducerReal.cxx:1381
lambda
Double_t lambda(Double_t x, Double_t y, Double_t z)
Definition: drawdal.C:48
PndDrcBarPoint::GetBarId
Int_t GetBarId() const
Definition: PndDrcBarPoint.h:43
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndGeoDrc::Lside
Double_t Lside()
Definition: PndGeoDrc.h:184
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
PndDrcHitProducerReal::SetPhotonDetEffNew
void SetPhotonDetEffNew(Double_t &lambda_min, Double_t &lambda_max, Double_t &lambda_step, Double_t efficiency[])
Definition: PndDrcHitProducerReal.cxx:691
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
PndGeoDrc::BarWidth
Double_t BarWidth()
Definition: PndGeoDrc.h:100
tr
Double_t tr
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:19
PndDrcHitProducerReal::FindDrcHitPosition
void FindDrcHitPosition(Double_t xPoint, Double_t yPoint, Double_t &xHit, Double_t &yHit, Int_t &pmtID)
Definition: PndDrcHitProducerReal.cxx:1345
n
int n
Definition: HypStatDecay.cxx:21
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
sigma
Double_t sigma[nsteps]
Definition: dedx_bands.C:65
PndDrcHitProducerReal::Init
virtual InitStatus Init()
Definition: PndDrcHitProducerReal.cxx:87
PndGeoDrc::BBoxGap
Double_t BBoxGap()
Definition: PndGeoDrc.h:130
PndGeoDrc::barBoxZDown
Double_t barBoxZDown()
Definition: PndGeoDrc.h:104
PndDrcBarPoint::GetThetaC
Double_t GetThetaC() const
Definition: PndDrcBarPoint.h:42
PndDrcHitProducerReal::SetPhotonTransportEff
void SetPhotonTransportEff(Double_t &, Double_t &, Double_t &, Double_t &, Int_t &, Double_t fEfficiency[])
Definition: PndDrcHitProducerReal.cxx:1301
PndDrcHitProducerReal::SetPhotonDetEffOld
void SetPhotonDetEffOld(Double_t &lambda_min, Double_t &lambda_max, Double_t &lambda_step, Double_t efficiency[])
Definition: PndDrcHitProducerReal.cxx:1239
pt
TString pt(TString pts, TString exts="px py pz")
Definition: invexp.C:133
PndGeoDrc::PipehAngle
Double_t PipehAngle()
Definition: PndGeoDrc.h:139
a
Int_t a
Definition: anaLmdDigi.C:126
PndGeoDrc::barBoxZUp
Double_t barBoxZUp()
Definition: PndGeoDrc.h:108
PndDrcBarPoint::GetBoxId
Int_t GetBoxId() const
Definition: PndDrcBarPoint.h:44
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
point
TClonesArray * point
Definition: anaLmdDigi.C:29
PndDrcHitProducerReal::SetFakeDetEff
void SetFakeDetEff(Double_t &lambda_min, Double_t &lambda_max, Double_t &lambda_step, Double_t efficiency[])
Definition: PndDrcHitProducerReal.cxx:668
PndDrcHitProducerReal::NumberOfBounces
void NumberOfBounces(TVector3, TVector3, Int_t *, Int_t *, Double_t *, Double_t *)
Definition: PndDrcHitProducerReal.cxx:439
PndDrcHitProducerReal::fGeo
PndGeoDrc * fGeo
Basic geometry data of barrel DRC.
Definition: PndDrcHitProducerReal.h:168
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndDrcHitProducerReal::AddPDHit
PndDrcPDHit * AddPDHit(Int_t detID, TVector3 posPDHit, TVector3 dPosPDHit, Double_t time, Double_t timeThreshold, Int_t index1)
Definition: PndDrcHitProducerReal.cxx:588
PndDrcHitProducerReal::Exec
virtual void Exec(Option_t *option)
Definition: PndDrcHitProducerReal.cxx:136
PndDrcHitProducerReal::FindPhiRot
Double_t FindPhiRot(Double_t, Double_t)
Definition: PndDrcHitProducerReal.cxx:515
hit
int hit(Int_t nEvents=0, TString inFile="sim.root", TString parFile="par.root", TString inDigi="digi.root", TString outFile="hit.root", Int_t timeBased=0)
Definition: hit.C:1
PndGeoDrc::barHalfThick
Double_t barHalfThick()
Definition: PndGeoDrc.h:96
P
GeV c P
Definition: ReadMesDecCorrRes_PiT.C:42
ClassImp
ClassImp(PndAnaContFact)
PndDrcHitProducerReal::AddHit
PndDrcHit * AddHit(Int_t detID, TVector3 posHit, TVector3 dPosHit, Double_t thetaC, Double_t errThetaC, Int_t index)
Definition: PndDrcHitProducerReal.cxx:572
t
TTree * t
Definition: bump_analys.C:13
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndDrcHitProducerReal::FindOutPoint
Double_t FindOutPoint(Double_t, Double_t, Double_t, Double_t *, Bool_t)
Definition: PndDrcHitProducerReal.cxx:538
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
NN
#define NN
Definition: TMVAClassification_BDT.class.C:81
PndGeoDrc::barNum
Double_t barNum()
Definition: PndGeoDrc.h:124
PndDrcHitProducerReal::FindDrcHitPositionTilt
void FindDrcHitPositionTilt(Double_t xPoint, Double_t yPoint, Double_t &xHit, Double_t &yHit, Int_t &pmtID)
Definition: PndDrcHitProducerReal.cxx:1362
thetaC
Double_t thetaC
Definition: plot_dirc.C:16
CAMath::Nint
int Nint(float x)
Definition: PndCAMath.h:117
PndGeoDrc::radius
Double_t radius()
Definition: PndGeoDrc.h:92
energy
Double_t energy
Definition: plot_dirc.C:15
Generated on Sun Aug 18 2019 03:01:20 for FairRoot/PandaRoot by   doxygen 1.8.8