FairRoot/PandaRoot
PndPhoGunShortP.cxx
Go to the documentation of this file.
1 
2 // -------------------------------------------------------------------------
3 // ----- PndPhoGunShortP source file -----
4 // ----- Created 11/09/13 by Harphool Kumawat -----
5 // ----- -----
6 // ----- -----
7 // -------------------------------------------------------------------------
8 #include <fstream>
9 #include <iostream>
10 #include "stdio.h"
11 
12 #include "PndGeoDrc.h"
13 #include "PndPhoGunShortP.h"
14 #include "FairRootManager.h"
15 #include "PndMCTrack.h"
16 #include "TVector3.h"
17 #include "PndDrcBarPoint.h"
18 #include "FairRunAna.h"
19 #include "FairRuntimeDb.h"
20 #include "FairBaseParSet.h"
21 #include "FairGeoVolume.h"
22 #include "TString.h"
23 #include "FairGeoTransform.h"
24 #include "FairGeoVector.h"
25 #include "FairGeoMedium.h"
26 #include "FairGeoNode.h"
27 #include "PndGeoDrcPar.h"
28 #include "TMath.h"
29 #include "TParticlePDG.h"
30 #include "TDatabasePDG.h"
31 #include "TPDGCode.h"
32 #include "TGeoManager.h"
33 #include "TFile.h"
34 #include "TExec.h"
35 
36 using std::endl;
37 using std::cout;
38 
39 // ----- Default constructor -------------------------------------------
40 PndPhoGunShortP::PndPhoGunShortP()
41 :FairTask("PndPhoGunShortP")
42 {
43  fGeo = new PndGeoDrc();
44  fGeoH=NULL;
45 // fTableName = "PndPhoGunShortP_table1.par";
46 // fOutputName = "PndPhoGunShortP_output1.root";
47 
48 }
49 // ----- Standard constructor with verbosity level -------------------------------------------
50 
51 PndPhoGunShortP::PndPhoGunShortP(Int_t verbose)
52  :FairTask("PndPhoGunShortP")
53 {
54  fVerbose = verbose;
55  fGeo = new PndGeoDrc();
56  fGeoH=NULL;
57 // fTableName = "PndPhoGunShortP_table1.par";
58 // fOutputName = "PndPhoGunShortP_output1.root";
59 }
60 // ----- Destructor ----------------------------------------------------
61 PndPhoGunShortP::~PndPhoGunShortP()
62 {
63  if (fGeo) delete fGeo;
64  if (fGeoH) delete fGeoH;
65  fHistoList->Delete();
66  delete fHistoList;
67 }
68 
69 // ----- Initialization -----------------------------------------------
70 InitStatus PndPhoGunShortP::Init()
71 {
72  cout << " ---------- INITIALIZATION ------------" << endl;
73  nevents = 0;
74  // Get RootManager
75  FairRootManager* ioman = FairRootManager::Instance();
76  if ( ! ioman ) {
77  cout << "-E- PndPhoGunShortP::Init: "
78  << "RootManager not instantiated!" << endl;
79  return kFATAL;
80  }
81 
82  // Get input array
83  fMCArray = (TClonesArray*) ioman->GetObject("MCTrack");
84  if ( ! fMCArray ) {
85  cout << "-W- PndPhoGunShortP::Init: "
86  << "No MCTrack array!" << endl;
87  return kERROR;
88  }
89  // Get Photon point array
90  fPDPointArray = (TClonesArray*) ioman->GetObject("DrcPDPoint");
91  if ( ! fPDPointArray ) {
92  cout << "-W- PndPhoGunShortP::Init: "
93  << "No DrcPDPoint array!" << endl;
94  return kERROR;
95  }
96  // Get input array
97  fPDHitArray = (TClonesArray*) ioman->GetObject("DrcPDHit");
98  if ( ! fPDHitArray ) {
99  cout << "-W- PndPhoGunShortP::Init: "
100  << "No DrcPDHit array!" << endl;
101  return kERROR;
102  }
103  fBarPointArray = (TClonesArray*) ioman->GetObject("DrcBarPoint");
104  if ( ! fBarPointArray ) {
105  cout << "-W- PndDrcLogLikeli::Init: "
106  << "No DrcBarPoint array!" << endl;
107  return kERROR;
108  }
109 
110  fEVPointArray = (TClonesArray*) ioman->GetObject("DrcEVPoint");
111  if ( ! fEVPointArray ) {
112  cout << "-W- PndPhoGunShortP::Init: "
113  << "No DrcEVPoint array!" << endl;
114  return kERROR;
115  }
116  fDigiArray = (TClonesArray*) ioman->GetObject("DrcDigi");
117  if ( ! fDigiArray ) {
118  cout << "-W- PndPhoGunShortP::Init: " << "No DrcDigi array!" << endl;
119  return kERROR;
120  }
121 
122  // Get parameters:
123  fpi = TMath::Pi();
124  fR = fGeo->radius();
125  fHThick = fGeo->barHalfThick();
126  fBboxNum = fGeo->BBoxNum();
127  fPipehAngle = fGeo->PipehAngle();
128  fBarBoxGap = fGeo->BBoxGap();
129  fLength = (180. - 2.*fPipehAngle - fBarBoxGap/fR*(fBboxNum/2. - 1.)/fpi*180.)/(fBboxNum/2.) * fR/ 180.*fpi;
130  fDphi = fGeo->BBoxAngle();//21.65 degrees = 2.*(180. - 2*fPipehAngle)/ fBboxNum; //[degrees]
131 
132  //Double_t EVlen = gGeoManager->GetVolume("DrcEVSensor")->GetShape()->Dz();
133  fEVlen = fGeo->EVlen();
134 
135  fEVdrop = fGeo->EVdrop();
136  fRBottom = fR - fHThick - fEVdrop + fR*(1.-cos(fDphi/2./180.*fpi))/cos(fDphi/2./180.*fpi);
137  //cout<<"fR = "<<fR<<", fHThick = "<<fHThick<<", fEVdrop = "<<fEVdrop<<", fRBottom = "<<fRBottom<<", dphi = "<<fDphi<<endl;
138 
139  Rin1 = fR - fHThick - fEVdrop - fGeo->boxGap() - fGeo->boxThick();
140  Rin2 = fR + fHThick + fGeo->EVoffset()+ fGeo->boxGap() + fGeo->boxThick() ;
141 // Rout1 = fR + fHThick + fGeo->EVoffset()+ fGeo->boxGap() + fGeo->boxThick() + 30.*tan(30.*fpi/180.);
142  Rout1 = fR - fHThick - fEVdrop -fGeo->boxGap() -fGeo->boxThick() + 5.*(fGeo->McpActiveArea() + fGeo->McpSize() - fGeo->McpActiveArea() + fGeo->McpGap()) -(fGeo->McpSize() - fGeo->McpActiveArea() + fGeo->McpGap());
143 
144  Ang_pipe = 180. - 2.*fPipehAngle;
145 
146  fFile = TFile::Open(fOutputName,"RECREATE");
147 // fLut = new TClonesArray("PndDrcLutNode");
148  fTree = new TTree("dircsim","Look-up table for DIRC");
149 // fTree->Branch("LUT",&fLut,256000,0);
150  for(Int_t iLut=0; iLut<5; iLut++){
151  fLut[iLut] = new TClonesArray("PndDrcLutNode");
152  fTree->Branch(Form("LUT%d",iLut),&fLut[iLut],256000,0);
153  }
154 
155  InitLut();
156 
157  if ( fGeoH == NULL ) fGeoH = PndGeoHandling::Instance();
158  fGeoH->SetParContainers();
159 
160  //fGeoH->GetGeoManager();
161  //const Double_t *truuu = gGeoManager->GetVolume("DrcEVSensor")->GetShape()->GetTransform()->GetTranslation();
162  //fLowZ = truuu[2]-fGeo->EVlen();
163  //cout<<"posZ = "<<fLowZ<<endl;
164 
165  cout << "-I- PndPhoGunShortP: Intialization successfull" << endl;
166  return kSUCCESS;
167 }
168 
169 // ----- Execution of Task ---------------------------------------------
170 void PndPhoGunShortP::Exec(Option_t*ion)
171 {
172  nevents++;
173 
174  fGeoH->SetVerbose(fVerbose);
175 
176  fDetectorID = 0;
177  if(nevents%1000==0) cout<<"Event # "<<nevents<<endl;
178  ProcessPhotonHit();
179 }
180 
181 //--------------Process Photon Hits----------------------------------------------------
182 void PndPhoGunShortP::ProcessPhotonHit()
183 {
184  Int_t EVEntry;
185  Int_t SelectionName=0;
186  Int_t PhiSec=0;
187  fmatrixdata.Set(9);
188  fmatrixdata.Reset(0);
189 
190  // Loop over PndDrcPDHits
191  for(Int_t k=0; k<fPDHitArray->GetEntriesFast(); k++) {
192 
193  pdhit = (PndDrcPDHit*)fPDHitArray->At(k);
194 
195  Int_t mcPDRef= pdhit->GetRefIndex();
196  fDigi = (PndDrcDigi*) fDigiArray->At(mcPDRef);
197  Int_t pointID= fDigi->GetIndex(0);
198  Ppt = (PndDrcPDPoint*)fPDPointArray->At(pointID);
199 
200 // Ppt = (PndDrcPDPoint*)fPDPointArray->At(mcPDRef);
201 
202  Int_t trID= Ppt->GetTrackID();
203  tr = (PndMCTrack*)fMCArray->At(trID);
204 
205  EVpt = (PndDrcEVPoint*)fEVPointArray->At(mcPDRef);
206 
207  if(trID!=EVpt->GetTrackID())cout<<"different track IDs"<<endl;
208  EVEntry = fEVPointArray->GetEntriesFast();
209  if(trID!=EVpt->GetTrackID())continue;
210 
211 // if(Ppt->GetBarPointID() < 0) continue;
212  PndDrcBarPoint *fBarPoint= (PndDrcBarPoint*)fBarPointArray->At(Ppt->GetBarPointID());
213 // fBarId = fBarPoint->GetDetectorID();
214  fBarId = fBarPoint->GetBarId();
215  TVector3 PdPoint;
216  PdPoint.SetXYZ(pdhit->GetX(),pdhit->GetY(),pdhit->GetZ());
217  Double_t PdPhi=InBarCoordSyst(PdPoint);
218 // if(fBarId<=0)cout<<"bar name - "<<fBarId<<endl;
219 // if(fBarId <= 0)continue;
220 // cout<<"boxname - "<<fBarPoint->GetBoxId()<<endl;
221  // production point of the photon
222  fStartVertex = tr->GetStartVertex();
223 
224  // find PHIrot to get to the bar coord. syst:
225  fPhiRot = InBarCoordSyst(fStartVertex);
226  fZin = ((PndDrcEVPoint*)fEVPointArray->At(0))->GetZ();
227  if(TMath::Nint(fabs(PdPhi-fPhiRot)*180./fpi/fDphi)>0)continue;
228 
229  // check if the last of the EV points is on the PD plane as the reflection was in the grease layer/window/photocathode, then exclude this point
230  if(((PndDrcEVPoint*)fEVPointArray->At(EVEntry-1))->GetZ() < fZin-fEVlen+0.01){EVEntry = EVEntry - 1;}
231  //check if there are EVpoints on the PD plane and these points are not the last in the array, then reject the photon
232  Bool_t wePho = kFALSE;
233  for(Int_t k2 = 1; k2 < fEVPointArray->GetEntriesFast()-1; k2++){
234  if((((PndDrcEVPoint*)fEVPointArray->At(k2))->GetZ()) < fZin-fEVlen+0.01){
235 // cout<<"this is a weird photon - has EV points at the PD plane!!!"<<endl;
236  fNweirdPhotons += 1;
237  wePho = kTRUE;
238  }
239  }
240  if(wePho == kTRUE)continue;
241 
242  if(EVEntry>5){
243  fNweirdPhotons += 1;
244  continue;
245  }
246 
247  if(EVEntry==1)ReflName="DD";
248  if(EVEntry>1){
249  ReflName="";
250  for(Int_t k1 = 1; k1 < EVEntry; k1++){
251  Int_t EVS=0;
252  EVt = (PndDrcEVPoint*)fEVPointArray->At(k1);
253  fEVSec.SetXYZ(EVt->GetX(),EVt->GetY(),EVt->GetZ());
254  fPhiRotEV=InBarCoordSyst(fEVSec);
255  fPhiRotEV=fPhiRotEV*180./fpi;
256  fEVPhi=fEVSec.Phi()*180./fpi;
257 // cout<<"EV phi = "<<fEVPhi<<", fPhiRot = "<<fPhiRot*180./fpi<<", fDphi = "<<fDphi<<endl;
258  PhiSec = TMath::Nint(fabs(fEVPhi-fPhiRot*180./fpi)/fDphi);
259  if(fEVPhi<0.){
260  fEVPhi=360.+fEVPhi;
261  PhiSec = TMath::Nint(fabs(fEVPhi-fPhiRot*180./fpi-360.)/fDphi);
262  }
263  if(PhiSec>0)cout<<"phi sec = "<<PhiSec<<endl;
264  FindReflectionType (EVt->GetX(),EVt->GetY(),EVt->GetZ(),ReflectionType);
265  ReflName.Append(ReflectionType);
266  }
267  if(PhiSec>0)continue;
268  if(EVEntry==2)ReflName.Append("1");
269  if(EVEntry==3)ReflName.Append("2");
270  if(EVEntry==4)ReflName.Append("3");
271  if(EVEntry==5)ReflName.Append("4");
272  if(EVEntry>=6)ReflName.Append("XX");
273 // ReflName.Prepend(Form("%d",PhiSec));
274  }
275  // 24 ambiguities in total:
276  if(ReflName == "DD"){fNoDD +=1; SelectionName=1;ambiguity =1;}
277  if(ReflName == "B1"){fNoB +=1; SelectionName=1;ambiguity =2;}
278  if(ReflName == "L1"){fNoL +=1; SelectionName=1;ambiguity =3;}
279  if(ReflName == "R1"){fNoR +=1; SelectionName=1;ambiguity =4;}
280  if(ReflName == "U1"){fNoU +=1; SelectionName=1;ambiguity =5;}
281  if(ReflName == "BU2"){fNoBU +=1; SelectionName=1;ambiguity =6;}
282  if(ReflName == "BL2"){fNoBL +=1; SelectionName=1;ambiguity =7;}
283  if(ReflName == "BR2"){fNoBR +=1; SelectionName=1;ambiguity =8;}
284  if(ReflName == "UL2"){fNoUL +=1; SelectionName=1;ambiguity =9;}
285  if(ReflName == "LU2"){fNoLU +=1; SelectionName=1;ambiguity =10;}
286  if(ReflName == "UR2"){fNoUR +=1; SelectionName=1;ambiguity =11;}
287  if(ReflName == "RU2"){fNoRU +=1; SelectionName=1;ambiguity =12;}
288  if(ReflName == "LR2"){fNoLR +=1; SelectionName=1;ambiguity =13;}
289  if(ReflName == "LB2"){fNoLB +=1; SelectionName=1;ambiguity =14;}
290  if(ReflName == "RB2"){fNoRB +=1; SelectionName=1;ambiguity =15;}
291  if(ReflName == "RL2"){fNoRL +=1; SelectionName=1;ambiguity =16;}
292  if(ReflName == "BLR3"){fNoBLR +=1; SelectionName=1;ambiguity =17;}
293  if(ReflName == "BUR3"){fNoBUR +=1; SelectionName=1;ambiguity =18;}
294  if(ReflName == "BUL3"){fNoBUL +=1; SelectionName=1;ambiguity =19;}
295  if(ReflName == "BLU3"){fNoBLU +=1; SelectionName=1;ambiguity =20;}
296  if(ReflName == "BRL3"){fNoBRL +=1; SelectionName=1;ambiguity =21;}
297  if(ReflName == "BRU3"){fNoBRU +=1; SelectionName=1;ambiguity =22;}
298  if(ReflName == "URL3"){fNoURL +=1; SelectionName=1;ambiguity =23;}
299  if(ReflName == "ULR3"){fNoULR +=1; SelectionName=1;ambiguity =24;}
300  if(ReflName == "RBL3"){fNoRBL +=1; SelectionName=1;ambiguity =25;}
301 // if(ReflName == "RBU3"){fNoRBU +=1; SelectionName=1;ambiguity =30;}
302  if(ReflName == "RLU3"){fNoRLU +=1; SelectionName=1;ambiguity =26;}
303 // if(ReflName == "RLB3"){fNoRLB +=1; SelectionName=1;ambiguity =32;}
304  if(ReflName == "RUL3"){fNoRUL +=1; SelectionName=1;ambiguity =27;}
305 // if(ReflName == "LRB3"){fNoLRB +=1; SelectionName=1;ambiguity =34;}
306  if(ReflName == "LRU3"){fNoLRU +=1; SelectionName=1;ambiguity =28;}
307  if(ReflName == "LBR3"){fNoLBR +=1; SelectionName=1;ambiguity =29;}
308  if(ReflName == "LBU3"){fNoLBU +=1; SelectionName=1;ambiguity =30;}
309  if(ReflName == "LUR3"){fNoLUR +=1; SelectionName=1;ambiguity =31;}
310 // if(ReflName == "LUB3"){fNoLUB +=1; SelectionName=1;ambiguity =39;}
311  if(ReflName == "URB3"){fNoURB +=1; SelectionName=1;ambiguity =32;}
312  if(ReflName == "ULB3"){fNoULB +=1; SelectionName=1;ambiguity =33;}
313 // if(ReflName == "UBR3"){fNoUBR +=1; SelectionName=1;ambiguity =42;}
314 // if(ReflName == "UBL3"){fNoUBL +=1; SelectionName=1;ambiguity =43;}
315  if(ReflName == "BULR4"){fNoBULR +=1; SelectionName=1;ambiguity =34;}
316  if(ReflName == "BURL4"){fNoBURL +=1; SelectionName=1;ambiguity =35;}
317  if(ReflName == "BLUR4"){fNoBLUR +=1; SelectionName=1;ambiguity =36;}
318  if(ReflName == "BRUL4"){fNoBRUL +=1; SelectionName=1;ambiguity =37;}
319  if(ReflName == "LBRU4"){fNoLBRU +=1; SelectionName=1;ambiguity =38;}
320  if(ReflName == "LBUR4"){fNoLBUR +=1; SelectionName=1;ambiguity =39;}
321  if(ReflName == "RBUL4"){fNoRBUL +=1; SelectionName=1;ambiguity =40;}
322  if(ReflName == "BLRU4"){fNoBLRU +=1; SelectionName=1;ambiguity =41;}
323  if(ReflName == "BRLU4"){fNoBRLU +=1; SelectionName=1;ambiguity =42;}
324  if(ReflName == "BURBL5"){fNoBURBL +=1; SelectionName=1;ambiguity =43;}
325  if(ReflName == "BRULB5"){fNoBRULB +=1; SelectionName=1;ambiguity =44;}
326  if(ReflName == "BULBR5"){fNoBULBR +=1; SelectionName=1;ambiguity =45;}
327  if(ReflName == "BURLB5"){fNoBURLB +=1; SelectionName=1;ambiguity =46;}
328  if(ReflName == "BLURB5"){fNoBLURB +=1; SelectionName=1;ambiguity =47;}
329  if(ReflName == "BULRB5"){fNoBULRB +=1; SelectionName=1;ambiguity =48;}
330 // if(ReflName == "BULB4"){fNoBULB +=1; SelectionName=1;ambiguity =46;}
331 // if(ReflName == "BURB4"){fNoBURB +=1; SelectionName=1;ambiguity =47;}
332 // if(ReflName == "UBLR4"){fNoUBLR +=1; SelectionName=1;ambiguity =48;}
333 // if(ReflName == "UBRL4"){fNoUBRL +=1; SelectionName=1;ambiguity =49;}
334 // if(ReflName == "ULBR4"){fNoULBR +=1; SelectionName=1;ambiguity =50;}
335 // if(ReflName == "URBL4"){fNoURBL +=1; SelectionName=1;ambiguity =51;}
336 // if(ReflName == "ULRB4"){fNoULRB +=1; SelectionName=1;ambiguity =52;}
337 // if(ReflName == "URLB4"){fNoURLB +=1; SelectionName=1;ambiguity =53;}
338 // if(ReflName == "RULR4"){fNoRULR +=1; SelectionName=1;ambiguity =54;}
339 // if(ReflName == "UBRL4"){fNoUBRL +=1; SelectionName=1;ambiguity =55;}
340 // if(ReflName == "RULB4"){fNoRULB +=1; SelectionName=1;ambiguity =56;}
341 // if(ReflName == "RBUL4"){fNoRBUL +=1; SelectionName=1;ambiguity =57;}
342 // if(ReflName == "RBLU4"){fNoRBLU +=1; SelectionName=1;ambiguity =58;}
343 // if(ReflName == "LURB4"){fNoLURB +=1; SelectionName=1;ambiguity =61;}
344 // if(ReflName == "LUBR4"){fNoLUBR +=1; SelectionName=1;ambiguity =62;}
345 // if(ReflName == "LRUL4"){fNoLRUL +=1; SelectionName=1;ambiguity =63;}
346 // if(ReflName == "LRLU4"){fNoLRLU +=1; SelectionName=1;ambiguity =64;}
347 // if(ReflName == "LURL4"){fNoLURL +=1; SelectionName=1;ambiguity =65;}
348 // if(ReflName == "RULR4"){fNoRULR +=1; SelectionName=1;ambiguity =66;}
349 // if(ReflName == "RBULR5"){fNoRBULR +=1; SelectionName=1;ambiguity =67;}
350 // if(ReflName == "RBLUR5"){fNoRBLUR +=1; SelectionName=1;ambiguity =68;}
351 // if(ReflName == "RBULB5"){fNoRBULB +=1; SelectionName=1;ambiguity =69;}
352 // if(ReflName == "RBRUL5"){fNoRBRUL +=1; SelectionName=1;ambiguity =70;}
353 // if(ReflName == "BRUBL5"){fNoBRUBL +=1; SelectionName=1;ambiguity =71;}
354 // if(ReflName == "BLRUL5"){fNoBLRUL +=1; SelectionName=1;ambiguity =74;}
355 // if(ReflName == "BLUBR5"){fNoBLUBR +=1; SelectionName=1;ambiguity =75;}
356 // if(ReflName == "LBURL5"){fNoLBURL +=1; SelectionName=1;ambiguity =80;}
357 // if(ReflName == "LBLUR5"){fNoLBLUR +=1; SelectionName=1;ambiguity =81;}
358 // cout<<"reflection type = "<<ReflName<<endl;
359  if(SelectionName==0.){continue;}
360 
361  fPixIndex = fDigi->GetSensorId();
362 // fPixIndex = pdhit->GetSensorId();
363 // fPixIndex = pdhit->GetDetectorID();
364  ftime = Ppt->GetTime();
365  // Photon initial momentum in the bar coord. syst.
366  fPphoInit = tr->GetMomentum();
367  //fPphoB = (fGeoH->MasterToLocalShortId(fPphoInit, fBarId) - fGeoH->MasterToLocalShortId((0.,0.,0.),fBarId)).Unit();
368  // fPphoB.Print();
369  // Photon initial momentum in the bar coord. syst.
370  fPphoB = (tr->GetMomentum()).Unit();
371  fPphoB.RotateZ(-fPhiRot);
372 // fkxBar = -fPphoB.Y();
373 // fkyBar = fPphoB.X();
374 // fkzBar = fPphoB.Z();
375  fkxBar = fPphoB.X();
376  fkyBar = fPphoB.Y();
377  fkzBar = fPphoB.Z();
378  fPphoB.SetXYZ(fkxBar,fkyBar,fkzBar);
379 // cout<<"pix "<<fPixIndex<<", reflection type = "<<ReflName<<", ambiguity = "<<ambiguity<<" barid "<<fBarId<<endl;
380  if(fPixIndex>30000)continue;
381 ((PndDrcLutNode*)(fLut[fBarId]->At(fPixIndex)))->AddEntry(fDigi->GetDetectorId(), fPphoB,ambiguity,ftime,pdhit->GetPosition());
382 
383 // ((PndDrcLutNodeH*)(fLut->At(fPixIndex)))->SetBarID(fBarId);
384 // ((PndDrcLutNodeH*)(fLut->At(fPixIndex)))->AddEntry(fPphoB,ambiguity,ftime);
385 // ((PndDrcLutNodeH*)(fLut->At(fPixIndex)))->SetPos(pdhit->GetPosition());
386 
387  fNoTotal +=1;
388 
389  }// photon hits
390 }
391 
392 //----------------------------------------------------------------------------------------------
393 Double_t PndPhoGunShortP::FindReflectionType(Double_t xev, Double_t yev, Double_t zev, TString ReflType){
394 
395 // cout<<"fPhiRotEV = "<<fPhiRotEV<< ", fDphi = "<<fDphi<<", Ang_pipe = "<<Ang_pipe<<endl;
396 
397  PlanB[0]=(Rin1/cos((Ang_pipe/16.)*fpi/180.))*cos((fPhiRotEV-fDphi/2.)*fpi/180.);
398  PlanB[1]=(Rin1/cos((Ang_pipe/16.)*fpi/180.))*sin((fPhiRotEV-fDphi/2.)*fpi/180.);
399  PlanB[2]=fZin;//-119.6015;//fGeo->barBoxZUp()-0.30075;
400  PlanB[3]=(Rin1/cos((Ang_pipe/16.)*fpi/180.))*cos((fPhiRotEV-fDphi/2.)*fpi/180.);
401  PlanB[4]=(Rin1/cos((Ang_pipe/16.)*fpi/180.))*sin((fPhiRotEV-fDphi/2.)*fpi/180.);
402  PlanB[5]=fZin-fEVlen;//-149.6015;//fGeo->barBoxZUp()-fGeo->EVlen()-0.30075;
403  PlanB[6]=(Rin1/cos((Ang_pipe/16.)*fpi/180.))*cos((fPhiRotEV+fDphi/2.)*fpi/180.);
404  PlanB[7]=(Rin1/cos((Ang_pipe/16.)*fpi/180.))*sin((fPhiRotEV+fDphi/2.)*fpi/180.);
405  PlanB[8]=fZin-fEVlen;//-149.6015;//fGeo->barBoxZUp()-fGeo->EVlen()-0.30075;
406 
407  PlanU[0]=(Rin2/cos((Ang_pipe/16.)*fpi/180.))*cos((fPhiRotEV-fDphi/2.)*fpi/180.);
408  PlanU[1]=(Rin2/cos((Ang_pipe/16.)*fpi/180.))*sin((fPhiRotEV-fDphi/2.)*fpi/180.);
409  PlanU[2]=fZin;//-119.6015;//fGeo->barBoxZUp()-0.30075;
410  PlanU[3]=(Rout1/cos((Ang_pipe/16.)*fpi/180.))*cos((fPhiRotEV-fDphi/2.)*fpi/180.);
411  PlanU[4]=(Rout1/cos((Ang_pipe/16.)*fpi/180.))*sin((fPhiRotEV-fDphi/2.)*fpi/180.);
412  PlanU[5]=fZin-fEVlen;//-149.6015;//fGeo->barBoxZUp()-fGeo->EVlen()-0.30075;
413  PlanU[6]=(Rout1/cos((Ang_pipe/16.)*fpi/180.))*cos((fPhiRotEV+fDphi/2.)*fpi/180.);
414  PlanU[7]=(Rout1/cos((Ang_pipe/16.)*fpi/180.))*sin((fPhiRotEV+fDphi/2.)*fpi/180.);
415  PlanU[8]=fZin-fEVlen;//-149.6015;//fGeo->barBoxZUp()-fGeo->EVlen()-0.30075;
416 
417  PlanR[0]=(Rin2/cos((Ang_pipe/16.)*fpi/180.))*cos((fPhiRotEV+fDphi/2.)*fpi/180.);
418  PlanR[1]=(Rin2/cos((Ang_pipe/16.)*fpi/180.))*sin((fPhiRotEV+fDphi/2.)*fpi/180.);
419  PlanR[2]=fZin;
420  PlanR[3]=(Rin1/cos((Ang_pipe/16.)*fpi/180.))*cos((fPhiRotEV+fDphi/2.)*fpi/180.);
421  PlanR[4]=(Rin1/cos((Ang_pipe/16.)*fpi/180.))*sin((fPhiRotEV+fDphi/2.)*fpi/180.);
422  PlanR[5]=fZin-fEVlen;
423 //cout<<Rin1<<" "<<Rin2<<" "<<Rout1<<" "<<Rout1/cos((Ang_pipe/16.)*fpi/180.)<<endl;
424  ReflectionType="";
425  TMatrixD matrix1;
426  fmatrixdata[0]=xev-PlanB[0];
427  fmatrixdata[1]=yev-PlanB[1];
428  fmatrixdata[2]=zev-PlanB[2];
429  fmatrixdata[3]=PlanB[3]-PlanB[0];
430  fmatrixdata[4]=PlanB[4]-PlanB[1];
431  fmatrixdata[5]=PlanB[5]-PlanB[2];
432  fmatrixdata[6]=PlanB[6]-PlanB[0];
433  fmatrixdata[7]=PlanB[7]-PlanB[1];
434  fmatrixdata[8]=PlanB[8]-PlanB[2];
435  matrix1.Use(3,3,fmatrixdata.GetArray());
436  determint1 = matrix1.Determinant();
437 
438  if(determint1 > -0.1 && determint1 < 0.1)ReflectionType="B";
439  fmatrixdata.Reset(0);
440  fmatrixdata[0]=xev-PlanU[0];
441  fmatrixdata[1]=yev-PlanU[1];
442  fmatrixdata[2]=zev-PlanU[2];
443  fmatrixdata[3]=PlanU[3]-PlanU[0];
444  fmatrixdata[4]=PlanU[4]-PlanU[1];
445  fmatrixdata[5]=PlanU[5]-PlanU[2];
446  fmatrixdata[6]=PlanU[6]-PlanU[0];
447  fmatrixdata[7]=PlanU[7]-PlanU[1];
448  fmatrixdata[8]=PlanU[8]-PlanU[2];
449  matrix1.Use(3,3,fmatrixdata.GetArray());
450  determint2 = matrix1.Determinant();
451  if(determint2 > -0.1 && determint2 < 0.1)ReflectionType="U";
452 
453  fmatrixdata.Reset(0);
454  fmatrixdata[0]=xev-PlanB[0];
455  fmatrixdata[1]=yev-PlanB[1];
456  fmatrixdata[2]=zev-PlanB[2];
457  fmatrixdata[3]=PlanU[3]-PlanB[0];
458  fmatrixdata[4]=PlanU[4]-PlanB[1];
459  fmatrixdata[5]=PlanU[5]-PlanB[2];;
460  fmatrixdata[6]=PlanB[3]-PlanB[0];
461  fmatrixdata[7]=PlanB[4]-PlanB[1];
462  fmatrixdata[8]=PlanB[5]-PlanB[2];
463  matrix1.Use(3,3,fmatrixdata.GetArray());
464  determint3 = matrix1.Determinant();
465 
466  fmatrixdata.Reset(0);
467  fmatrixdata[0]=xev-PlanR[0];
468  fmatrixdata[1]=yev-PlanR[1];
469  fmatrixdata[2]=zev-PlanR[2];
470  fmatrixdata[3]=PlanR[3]-PlanR[0];
471  fmatrixdata[4]=PlanR[4]-PlanR[1];
472  fmatrixdata[5]=PlanR[5]-PlanR[2];
473  fmatrixdata[6]=PlanU[6]-PlanR[0];
474  fmatrixdata[7]=PlanU[7]-PlanR[1];
475  fmatrixdata[8]=PlanU[8]-PlanR[2];
476  matrix1.Use(3,3,fmatrixdata.GetArray());
477  determint4 = matrix1.Determinant();
478 
479  if(ReflectionType!="U" && ReflectionType!="B" && determint3 > determint4)ReflectionType="R";
480  if(ReflectionType!="U" && ReflectionType!="B" && determint3 < determint4)ReflectionType="L";
481 // cout<<"determinants: "<<determint1<<", "<<determint2<<endl;
482 // cout<<"determinants: "<<determint3<<", "<<determint4<<endl;
483 // cout<<"rin2 = "<<Rin1<<", rin2 = "<<Rin2<<", rout = "<<Rout1<<endl;
484 }
485 
486 //----------------------------------------------------------------------------------------------
487 Double_t PndPhoGunShortP::InBarCoordSyst(TVector3 start){
488 
489  // this function is used with fStartVertex to find the bar from which the photon originated
490 
491  Double_t startPhi = start.Phi()/fpi*180.; // [degrees]
492  //cout<<"-I- InBarCoordinateSystem: start phi = "<<startPhi<<endl;
493  //cout<<"-I- InBarCoordinateSystem: dphi = "<<Dphi<<endl;
494  Double_t PhiRot = 0.; //[degrees]
495  if(startPhi < 0.){startPhi = 360. + startPhi;}
496  if(startPhi > 0. && startPhi < 90.){
497  PhiRot = TMath::Floor(startPhi/fDphi) *fDphi + fDphi/2.;
498  }
499  if(startPhi > 90. && startPhi < 270.){
500  PhiRot = 90. + fPipehAngle + TMath::Floor((startPhi-90.-fPipehAngle)/fDphi) *fDphi + fDphi/2.;
501  }
502  if(startPhi > 270. && startPhi < 360.){
503  PhiRot = 270. + fPipehAngle + TMath::Floor((startPhi-270.-fPipehAngle)/fDphi) *fDphi + fDphi/2.;
504  }
505  //cout<<"-I- InBarCoordinateSystem: PhiRot = "<<PhiRot<<endl;
506 
507 
508  return PhiRot/180.*fpi;
509 }
510 //--------------------------------------------------------------------------------------------
511 void PndPhoGunShortP::InitLut()
512 
513 {
514  Int_t Nnodes = 30000;
515  for(Int_t iLut=0; iLut<5; iLut++){
516  TClonesArray &fLuta = *fLut[iLut];
517  for (Long64_t n=0; n<Nnodes; n++) {
518  new((fLuta)[n]) PndDrcLutNode(-1);
519  }
520  }
521 }
522 
523 // ----- Finish Task ---------------------------------------------------
524 void PndPhoGunShortP::Finish()
525 {
526  fTree->Fill();
527  fTree->Write();
528  fFile->Write();
529 
530  for(Int_t iLut=0; iLut<5; iLut++){
531  fLut[iLut]->Clear();
532  }
533  cout << "-I- PndDrcLutFill: Finish" << endl;
534 }
535 
536 // -------------------------------------------------------------------------
537 ClassImp(PndPhoGunShortP)
PndPhoGunShortP::pdhit
PndDrcPDHit * pdhit
Definition: PndPhoGunShortP.h:153
PndPhoGunShortP::fEVPointArray
TClonesArray * fEVPointArray
Definition: PndPhoGunShortP.h:82
PndPhoGunShortP::~PndPhoGunShortP
virtual ~PndPhoGunShortP()
Definition: PndPhoGunShortP.cxx:61
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndGeoDrc::BBoxNum
Double_t BBoxNum()
Definition: PndGeoDrc.h:136
PndDrcBarPoint::GetBarId
Int_t GetBarId() const
Definition: PndDrcBarPoint.h:43
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
PndPhoGunShortP::fDigiArray
TClonesArray * fDigiArray
Definition: PndPhoGunShortP.h:83
PndGeoHandling::SetParContainers
virtual void SetParContainers()
Definition: PndGeoHandling.cxx:72
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
PndGeoDrc::EVlen
Double_t EVlen()
Definition: PndGeoDrc.h:127
n
int n
Definition: HypStatDecay.cxx:21
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndPhoGunShortP::EVpt
PndDrcEVPoint * EVpt
Definition: PndPhoGunShortP.h:150
PndPhoGunShortP::fPDHitArray
TClonesArray * fPDHitArray
Definition: PndPhoGunShortP.h:80
PndPhoGunShortP::fGeo
PndGeoDrc * fGeo
Basic geometry data of barrel DRC.
Definition: PndPhoGunShortP.h:91
PndPhoGunShortP::InBarCoordSyst
Double_t InBarCoordSyst(TVector3)
Definition: PndPhoGunShortP.cxx:487
PndGeoDrc::BBoxGap
Double_t BBoxGap()
Definition: PndGeoDrc.h:130
PndGeoDrc::EVdrop
Double_t EVdrop()
Definition: PndGeoDrc.h:142
PndGeoDrc::boxGap
Double_t boxGap()
Definition: PndGeoDrc.h:116
PndGeoDrc::EVoffset
Double_t EVoffset()
Definition: PndGeoDrc.h:145
PndGeoDrc::PipehAngle
Double_t PipehAngle()
Definition: PndGeoDrc.h:139
PndDrcPDPoint::GetBarPointID
Int_t GetBarPointID() const
Definition: PndDrcPDPoint.h:56
PndDrcPDHit::GetRefIndex
virtual Int_t GetRefIndex()
Definition: PndDrcPDHit.h:57
PndPhoGunShortP::Finish
virtual void Finish()
Definition: PndPhoGunShortP.cxx:524
PndPhoGunShortP::fNweirdPhotons
Double_t fNweirdPhotons
Definition: PndPhoGunShortP.h:130
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndGeoDrc::McpActiveArea
Double_t McpActiveArea()
Definition: PndGeoDrc.h:166
PndGeoDrc::McpSize
Double_t McpSize()
Definition: PndGeoDrc.h:163
PndGeoDrc::BBoxAngle
Double_t BBoxAngle()
Definition: PndGeoDrc.h:133
PndPhoGunShortP::EVt
PndDrcEVPoint * EVt
Definition: PndPhoGunShortP.h:151
PndPhoGunShortP::fLut
TClonesArray * fLut[5]
Definition: PndPhoGunShortP.h:155
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndPhoGunShortP::fGeoH
PndGeoHandling * fGeoH
Definition: PndPhoGunShortP.h:92
PndGeoDrc::boxThick
Double_t boxThick()
Definition: PndGeoDrc.h:120
PndGeoHandling::Instance
static PndGeoHandling * Instance()
Definition: PndGeoHandling.cxx:36
PndDrcPDHit::GetPosition
TVector3 GetPosition() const
Definition: PndDrcPDHit.h:58
PndPhoGunShortP::FindReflectionType
Double_t FindReflectionType(Double_t, Double_t, Double_t, TString)
Definition: PndPhoGunShortP.cxx:393
PndPhoGunShortP::fPDPointArray
TClonesArray * fPDPointArray
Definition: PndPhoGunShortP.h:79
PndGeoHandling::SetVerbose
void SetVerbose(Int_t v)
Definition: PndGeoHandling.h:125
PndPhoGunShortP::fMCArray
TClonesArray * fMCArray
Definition: PndPhoGunShortP.h:81
PndPhoGunShortP::PlanR
Double_t PlanR[6]
Definition: PndPhoGunShortP.h:139
PndGeoDrc::barHalfThick
Double_t barHalfThick()
Definition: PndGeoDrc.h:96
PndPhoGunShortP::Init
virtual InitStatus Init()
Definition: PndPhoGunShortP.cxx:70
PndDrcDigi::GetSensorId
Int_t GetSensorId() const
Definition: PndDrcDigi.h:91
PndPhoGunShortP::Exec
virtual void Exec(Option_t *option)
Definition: PndPhoGunShortP.cxx:170
ClassImp
ClassImp(PndAnaContFact)
PndPhoGunShortP::Ppt
PndDrcPDPoint * Ppt
Definition: PndPhoGunShortP.h:152
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndPhoGunShortP::tr
PndMCTrack * tr
Definition: PndPhoGunShortP.h:154
PndPhoGunShortP::PlanU
Double_t PlanU[9]
Definition: PndPhoGunShortP.h:139
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndDrcDigi::GetDetectorId
Int_t GetDetectorId() const
Definition: PndDrcDigi.h:90
At
cout<<"the Event No is "<< i<< endl;{{if(hit_array->GetEntriesFast()!=mc_array->GetEntriesFast()) continue;PndSdsHit *hit=(PndSdsHit *) hit_array-> At(j)
Definition: anaLmdCluster.C:71
PndGeoDrc::McpGap
Double_t McpGap()
Definition: PndGeoDrc.h:169
CAMath::Nint
int Nint(float x)
Definition: PndCAMath.h:117
PndPhoGunShortP::PlanB
Double_t PlanB[9]
Definition: PndPhoGunShortP.h:139
PndDrcDigi::GetIndex
Int_t GetIndex(int i=0) const
Definition: PndDrcDigi.h:94
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
PndGeoDrc::radius
Double_t radius()
Definition: PndGeoDrc.h:92
PndPhoGunShortP::fDigi
PndDrcDigi * fDigi
Definition: PndPhoGunShortP.h:90
PndPhoGunShortP::fBarPointArray
TClonesArray * fBarPointArray
Definition: PndPhoGunShortP.h:78
Generated on Wed Apr 3 2019 10:02:55 for FairRoot/PandaRoot by   doxygen 1.8.5