#include <PndDrcRecoLookupMapS.h>

Inheritance diagram for PndDrcRecoLookupMapS:

Public Member Functions
	PndDrcRecoLookupMapS ()

	PndDrcRecoLookupMapS (Int_t verbose)

virtual	~PndDrcRecoLookupMapS ()

virtual InitStatus	Init ()

virtual void	Exec (Option_t *option)

virtual void	Finish ()

void	SetHistoWidth (double wid=0.01)

void	SetNSigma (double nsig=3.)

void	SetMagneticField (double mag=0.)

void	SetOutputFile (TString outName="PndDrcRecoLookupMapS_output.root")

Private Member Functions
void	ProcessPhotonHit ()

void	ProcessPhotonMC ()

void	WriteToFile ()

void	CreateHisto ()

void	DrawHisto ()

void	ProcessBarHit ()

virtual void	SetParContainers ()

Double_t	InBarCoordSyst (TVector3, TVector3 , TVector3 , TVector3 , TVector3 )

Int_t	NumberOfBounces (TVector3, TVector3, Int_t)

Double_t	FindOutPoint (Double_t, Double_t, Double_t, Double_t *, Bool_t)

Double_t	CherenkovPhi (TVector3, TVector3, Double_t , Double_t )

Double_t	RecoAmbigTime (TVector3, TVector3, Double_t *, Bool_t)

void	DrawBarBox (TVector3, TVector3, TVector3, TVector3)

void	SetDefaultParameters ()

Private Attributes
Double_t	fpi

Double_t	fR

Double_t	fzup

Double_t	fzdown

Double_t	fHThick

Double_t	fBboxNum

Double_t	fBarNum

Double_t	fPipehAngle

Double_t	fBarBoxGap

Double_t	fLSide

Double_t	fBarWidth

Double_t	fDphi

TTree *	mapstree

Float_t	fxx

Float_t	fyy

Float_t	fchdiff

PndDrcDigiPar *	fDigiPar

Int_t	fDetectorID

TClonesArray *	fBarPointArray

TClonesArray *	fPDPointArray

TClonesArray *	fHitArray

TClonesArray *	fPDHitArray

TClonesArray *	fMCArray

TClonesArray *	fDigiArray

TClonesArray *	fDrcLutInfoArray

TList *	fHistoList

PndGeoDrc *	fGeo
	Basic geometry data of barrel DRC. More...

PndGeoHandling *	fGeoH

Int_t	fVerbose

Int_t	nevents

TString	fOutputName

Double_t	fWidth

Double_t	fPixelSize

Double_t	fHAngleInBDeg

Double_t	fCHrecoMin

Double_t	fCHrecoMax

Double_t	fNSigma

Double_t	fB

TVector3	fnX1

TVector3	fnY1

Double_t	flambdah

Double_t	fpixID

Double_t	fxbar

Double_t	fybar

Double_t	fxPHit

Double_t	fyPHit

Double_t	fzPHit

Double_t	ftime

Double_t	ftime0

Double_t	fPhiRot

TVector3	fPMo

TVector3	fPMoB

TVector3	fPphoInit

TVector3	fkBar

Double_t	fWeight

Int_t	fBarId

Bool_t	print

Int_t	NHitPix

Int_t	NPixPar

Int_t	NAmb

Double_t	fkxBar

Double_t	fkyBar

Double_t	fkzBar

Double_t	fPxMoBar

Double_t	fPyMoBar

Double_t	fPzMoBar

Double_t	fxnew

Double_t	fynew

Double_t	CHexp

Double_t	CHreal

TVector3	fPphoB

TVector3	fStartVertex

TVector3	fPphoPD

Double_t	fPath

Double_t	fPhiMap

Double_t	fThetaMap

TVector3	fBBver1

TVector3	fBBver2

TVector3	fBBver3

TVector3	fBBver4

Detailed Description

PndDrcRecoLookupMapS.h

Class for Analysising DRC Cherenkov Photons with lookup tables for CHerenkov photons directions

Definition at line 49 of file PndDrcRecoLookupMapS.h.

Constructor & Destructor Documentation

PndDrcRecoLookupMapS::PndDrcRecoLookupMapS ( )

Default constructor

Definition at line 55 of file PndDrcRecoLookupMapS.cxx.

References fDigiPar, fGeo, and fGeoH.

 :FairTask("PndDrcRecoLookupMapS")
 {
   fGeo = new PndGeoDrc();
   fGeoH=NULL;
   fDigiPar = NULL;  
 }

PndDrcRecoLookupMapS::PndDrcRecoLookupMapS ( Int_t verbose )

Constructor with verbosity

Definition at line 64 of file PndDrcRecoLookupMapS.cxx.

References fDigiPar, fGeo, fGeoH, fVerbose, and verbose.

   :FairTask("PndDrcRecoLookupMapS")
 {
   fVerbose = verbose;  
   fGeo = new PndGeoDrc();
   fGeoH=NULL;
   fDigiPar = NULL;
 }

PndDrcRecoLookupMapS::~PndDrcRecoLookupMapS ( )

virtual

Destructor

Definition at line 73 of file PndDrcRecoLookupMapS.cxx.

References fGeo, fGeoH, and fHistoList.

 {
   if (fGeo) delete fGeo;
   if (fGeoH) delete fGeoH;
   fHistoList->Delete();  
   delete fHistoList;      
 }

Member Function Documentation

Double_t PndDrcRecoLookupMapS::CherenkovPhi	(	TVector3	,
		TVector3	,
		Double_t *	,
		Double_t *
	)

private

void PndDrcRecoLookupMapS::CreateHisto ( )

private

Definition at line 690 of file PndDrcRecoLookupMapS.cxx.

Referenced by Init().

690 {}

void PndDrcRecoLookupMapS::DrawBarBox	(	TVector3	v1,
		TVector3	v2,
		TVector3	v3,
		TVector3	v4
	)

private

Definition at line 633 of file PndDrcRecoLookupMapS.cxx.

References Double_t.

                                                                                        {
   Double_t xx[5];
   Double_t yy[5];
   xx[0] = v1.X(); xx[1] = v2.X(); xx[2] = v3.X(); xx[3] = v4.X(); xx[4] = v1.X();       
   yy[0] = v1.Y(); yy[1] = v2.Y(); yy[2] = v3.Y(); yy[3] = v4.Y(); yy[4] = v1.Y();
   TPolyLine *p6 =new TPolyLine(5,xx,yy);
   p6->SetLineWidth(2);
   p6->Draw("same");
 }

void PndDrcRecoLookupMapS::DrawHisto ( )

private

Definition at line 704 of file PndDrcRecoLookupMapS.cxx.

Referenced by Finish().

704 {}

void PndDrcRecoLookupMapS::Exec ( Option_t * option )

virtual

Executed task

Definition at line 196 of file PndDrcRecoLookupMapS.cxx.

References fDetectorID, fGeoH, fVerbose, nevents, ProcessPhotonHit(), and PndGeoHandling::SetVerbose().

 {
   //if ( ! fChPhoArray ) Fatal("Exec", "No fChPhoArray");
   //fChPhoArray->Clear();
   fGeoH->SetVerbose(fVerbose);
   
   nevents++;
   fDetectorID = 0;
    
   cout<<"EVENT # "<<nevents<<endl;
   // fNHits = 0;
   //ProcessBarHit();
   ProcessPhotonHit();
   //ProcessPhotonMC();
 }

Double_t PndDrcRecoLookupMapS::FindOutPoint	(	Double_t	x0,
		Double_t	xEn,
		Double_t	a,
		Double_t *	NN,
		Bool_t	printt
	)

private

Definition at line 599 of file PndDrcRecoLookupMapS.cxx.

References a, Double_t, m, and n.

Referenced by NumberOfBounces().

                                                                                                              {  
         Double_t m=99.;
         Double_t n=TMath::Floor(x0/a);
         m = n;          
         if(printt){std::cout<<"n = "<<n<<", NN = "<<*NN<<", x0 = "<<x0<<", a = "<<a<<std::endl;}
         Double_t x1 = x0 - n*a;
         if(x0 < 0.){x1 = x0 - (n+1)*a;}
         if(printt){std::cout<<"xy = "<< x1<<std::endl;}
         Double_t xK = 0.;
         if((m/2. - TMath::Floor(m/2.)) == 0.) { // 4etnoe
                 //if(print){std::cout<<"odd==0"<<std::endl;}
                 if(x0 >= 0. && x1 + xEn <= a){xK = x1 + xEn;} 
                 if(x0 >= 0. && x1 + xEn >  a){xK = 2*a - x1 - xEn; n = 1. + n;}
                 if(x0 < 0. && x1 + xEn >= 0.){xK = a - (x1 + xEn); n = -1. -n;}
                 if(x0 < 0. && x1 + xEn < 0.) {xK = a + x1 + xEn; n = -n;}
         if(printt){std::cout<<"xK = "<< xK<<", n = "<<n<<std::endl;}
                 
         }
 
         if((m/2. - TMath::Floor(m/2.)) != 0.) { // ne4etnoe
                 //if(print){std::cout<<"even!=0"<<std::endl;} 
                 if(x0 >= 0. && x1 + xEn <= a){xK = a - (x1 + xEn);} 
                 if(x0 >= 0. && x1 + xEn >  a){xK = x1 + xEn - a; n = 1. + n;}
                 if(x0 < 0. && x1 + xEn >= 0.){xK = x1 + xEn; n = -1. -n;} 
                 if(x0 < 0. && x1 + xEn < 0.) {xK = - (x1 + xEn); n = -n;}
         if(printt){std::cout<<"xK = "<< xK<<", n = "<<n<<std::endl;}
 
         }
         
         *NN = n;                
         return xK;
 }

void PndDrcRecoLookupMapS::Finish ( )

virtual

Finish task

Definition at line 696 of file PndDrcRecoLookupMapS.cxx.

References DrawHisto(), and WriteToFile().

 {
   cout << "-I- PndDrcRecoLookupMapS: Finish" << endl;  
   WriteToFile();
   DrawHisto();
 }

Double_t PndDrcRecoLookupMapS::InBarCoordSyst	(	TVector3	start,
		TVector3 *	v1,
		TVector3 *	v2,
		TVector3 *	v3,
		TVector3 *	v4
	)

private

Definition at line 523 of file PndDrcRecoLookupMapS.cxx.

References Double_t, fDphi, fHThick, fLSide, fpi, fPipehAngle, and fR.

Referenced by ProcessPhotonHit().

                                                                                                                    {
     Double_t startPhi = start.Phi()/fpi*180.; // [degrees]
     //cout<<"-I- InBarCoordinateSystem: start phi = "<<start.Phi()/fpi*180.<<endl;
     if(startPhi < 0.){startPhi = 360. + startPhi;}
     //if(startPhi > 360.){startPhi = startPhi - 360.;}
     //cout<<"-I- InBarCoordinateSystem: start phi = "<<startPhi<<endl;
     //Double_t dphi = 2.*(180. - 2*fPipehAngle)/ fBboxNum; //[degrees]
     //cout<<"-I- InBarCoordinateSystem: dphi = "<<dphi<<endl;
     Double_t PhiRot = 0.; //[degrees]
     if(startPhi >= 0. && startPhi < 90.){
       PhiRot = TMath::Floor(startPhi/fDphi) *fDphi + fDphi/2.;
     }
     if(startPhi >= 90. && startPhi < 270.){
       PhiRot = 90.  + fPipehAngle + TMath::Floor((startPhi-90.-fPipehAngle)/fDphi) *fDphi + fDphi/2.;
     } 
     if(startPhi >= 270. && startPhi < 360.){
       PhiRot = 270. + fPipehAngle + TMath::Floor((startPhi-270.-fPipehAngle)/fDphi) *fDphi + fDphi/2.;
     }
     //cout<<"-I- InBarCoordinateSystem: PhiRot = "<<PhiRot<<endl;
       
     TVector3 ver1, ver2, ver3, ver4;
     ver1.SetXYZ(fR-(fHThick), fLSide/2. ,0.);
     ver2.SetXYZ(fR+(fHThick), fLSide/2. ,0.);
     ver3.SetXYZ(fR+(fHThick), -fLSide/2. ,0.);
     ver4.SetXYZ(fR-(fHThick), -fLSide/2. ,0.);
     
     ver1.RotateZ(PhiRot/180.*fpi);
     ver2.RotateZ(PhiRot/180.*fpi);
     ver3.RotateZ(PhiRot/180.*fpi);
     ver4.RotateZ(PhiRot/180.*fpi);
     
     *v1 = ver1;
     *v2 = ver2;
     *v3 = ver3;
     *v4 = ver4;   
     
     return PhiRot/180.*fpi;
 }

InitStatus PndDrcRecoLookupMapS::Init ( )

virtual

Definition at line 82 of file PndDrcRecoLookupMapS.cxx.

References PndGeoDrc::barBoxZDown(), PndGeoDrc::barBoxZUp(), PndGeoDrc::barHalfThick(), PndGeoDrc::barNum(), PndGeoDrc::BarWidth(), PndGeoDrc::BBoxGap(), PndGeoDrc::BBoxNum(), CreateHisto(), fBarBoxGap, fBarNum, fBarPointArray, fBarWidth, fBboxNum, fDigiArray, fDphi, fDrcLutInfoArray, fGeo, fHThick, fLSide, fMCArray, fnX1, fnY1, fPDHitArray, fPDPointArray, fpi, fPipehAngle, fR, fzdown, fzup, PndGeoDrc::Lside(), nevents, Pi, PndGeoDrc::PipehAngle(), and PndGeoDrc::radius().

 {
   cout << " ---------- INITIALIZATION ------------" << endl;
   nevents = 0;
   // Get RootManager
   FairRootManager* ioman = FairRootManager::Instance();
   if ( ! ioman ) {
     cout << "-E- PndDrcRecoLookupMapS::Init: "
          << "RootManager not instantiated!" << endl;
     return kFATAL;
   }
 
   // Get input array
   fMCArray = (TClonesArray*) ioman->GetObject("MCTrack");
   if ( ! fMCArray ) {
     cout << "-W- PndDrcRecoLookupMapS::Init: "
          << "No MCTrack array!" << endl;
     return kERROR;
   }
   // Get input array
   fBarPointArray = (TClonesArray*) ioman->GetObject("DrcBarPoint");
   if ( ! fBarPointArray ) {
     cout << "-W- PndDrcRecoLookupMapS::Init: "
          << "No DrcBarPoint array!" << endl;
     return kERROR;
   }
   // Get Photon point array
   fPDPointArray = (TClonesArray*) ioman->GetObject("DrcPDPoint");
   if ( ! fPDPointArray ) {
     cout << "-W- PndDrcRecoLookupMapS::Init: "
          << "No DrcPDPoint array!" << endl;
     return kERROR;
   }
 /*  // Get input array
   fHitArray = (TClonesArray*) ioman->GetObject("DrcHit");
   if ( ! fHitArray ) {
     cout << "-W- PndDrcRecoLookupMapS::Init: "
          << "No DrcHit array!" << endl;
          return kERROR;
          }
 */
 
   // Get digi array
   fDigiArray = (TClonesArray*) ioman->GetObject("DrcDigi");
   if ( ! fDigiArray ) {
     cout << "-W- PndDrcRecoLookupMapS::Init: " << "No DrcDigi array!" << endl;
     return kERROR;
   }
 
   // Get input array
   fPDHitArray = (TClonesArray*) ioman->GetObject("DrcPDHit");
   if ( ! fPDHitArray ) {
     cout << "-W- PndDrcRecoLookupMapS::Init: "
          << "No DrcPDHit array!" << endl;
     return kERROR;
   }
     
    // Create and register output array
   // fChPhoArray = new TClonesArray("PndChPho");
   // ioman->Register("ChPho","Drc",fChPhoArray, kTRUE);  
   fDrcLutInfoArray = new TClonesArray("PndDrcLutInfo");
   ioman->Register("DrcLutInfo","Drc",fDrcLutInfoArray, kTRUE); 
   
    // Get parameters:
    fpi = TMath::Pi();
    
    fR          = fGeo->radius();
    fzup        = fGeo->barBoxZUp();
    fzdown      = fGeo->barBoxZDown();
    fHThick     = fGeo->barHalfThick();
    fBboxNum    = fGeo->BBoxNum();
    fBarNum     = fGeo->barNum();
    fPipehAngle = fGeo->PipehAngle();
    fBarBoxGap  = fGeo->BBoxGap();
    //fLength     = (180. - 2.*fPipehAngle - fBarBoxGap/fR*(fBboxNum/2. - 1.)/fpi*180.)/(fBboxNum/2.) * fR/ 180.*fpi;
    fDphi       = 2.*(180. - 2*fPipehAngle)/ fBboxNum; //[degrees]
    fLSide      = fGeo->Lside();//(180. - 2.*fPipehAngle - fBarBoxGap/fR*(fBboxNum/2. - 1.)/fpi*180.)/(fBboxNum/2.) * fR/ 180.*fpi;
    fBarWidth   = fGeo->BarWidth();//fLSide/fBarNum;
    
    // vectors of bar surfaces:
    fnX1.SetXYZ(-1., 0.,0.);   
    fnY1.SetXYZ( 0., 1.,0.);   
   
    cout << "-I- PndDrcRecoLookupMapS: Intialization successfull" << endl;
    CreateHisto();
    return kSUCCESS;
    
 }

Int_t PndDrcRecoLookupMapS::NumberOfBounces	(	TVector3	start,
		TVector3	dir,
		Int_t	barId
	)

private

Definition at line 563 of file PndDrcRecoLookupMapS.cxx.

References CAMath::Cos(), Double_t, fabs(), fBarWidth, fGeoH, fHThick, FindOutPoint(), fzdown, fzup, PndGeoHandling::MasterToLocalShortId(), CAMath::Sin(), and CAMath::Tan().

Referenced by ProcessPhotonHit().

                                                                                     {
     // start - photon production point in global coord system
     // dir - photon direction in bar coord system
     
     //cout<<"-I- NumberOfBounces: dir"<<endl;
     //dir.Print();
     //cout<<"-I- NumberOfBounces: start"<<endl;
     //start.Print();   
        
     // Find coordinates of X0, Y0:
     Double_t Z0, X0, Y0;
     if(dir.Theta() < 3.1415/2.){
       Z0 = -(fabs(fzup) + 2.*fzdown - start.Z());
     }
     if(dir.Theta() >= 3.1415/2.){
       Z0 = -(start.Z() -  fzup);
     }    
     X0 = Z0*TMath::Tan(dir.Theta())*TMath::Cos(dir.Phi());
     Y0 = Z0*TMath::Tan(dir.Theta())*TMath::Sin(dir.Phi());
     //cout<<"-I- NumberOfBounces: tan th = "<<TMath::Tan(dir.Theta())<<", sin ph = "<<TMath::Sin(dir.Phi())<<", cos ph = "<<TMath::Cos(dir.Phi())<<endl;
     //cout<<"-I- NumberOfBounces: X0 = "<<X0<<", Y0 = "<<Y0<<", Z0 = "<<Z0<<endl;
     
     // Find the start position of the photon with respect to the middle of the bar:
     TVector3 startLocal = fGeoH->MasterToLocalShortId(start, barId); // point
     //cout<<"-I- NumberOfBounces: Xen = "<<startLocal.X() + fBarWidth/2.<<", Yen = "<<startLocal.Y() + fHThick<<endl;
         
     // Find the number of bounces in each direction       
     Double_t N1, N2;   
     FindOutPoint(X0, startLocal.X() + fBarWidth/2., fBarWidth,  &N1, 0);
     FindOutPoint(Y0, startLocal.Y() + fHThick,      2.*fHThick, &N2, 0);
     //cout<<"-I- NumberOfBounces: N1 = "<<N1<<", N2 = "<<N2<<endl;
    
     return (Int_t)N1+N2;
 }

void PndDrcRecoLookupMapS::ProcessBarHit ( )

private

Definition at line 223 of file PndDrcRecoLookupMapS.cxx.

References fBarPointArray, fHitArray, fMCArray, PndDrcHit::GetRefIndex(), hit, pt(), and tr.

 {
  PndDrcBarPoint* pt=NULL;
   PndDrcHit* hit=NULL;
   PndMCTrack* tr = NULL;
   
   for(Int_t j=0; j<fHitArray->GetEntriesFast(); j++) {
     hit = (PndDrcHit*)fHitArray->At(j);
     Int_t mcRef= hit->GetRefIndex();
     pt= (PndDrcBarPoint*)fBarPointArray->At(mcRef);
     
     Int_t chtrID= pt->GetTrackID();
     tr = (PndMCTrack*)fMCArray->At(chtrID);
     
     cout<<"nother track px = "<<pt->GetPx()<<", py = "<<pt->GetPy()<<endl;
     
   }
 }

void PndDrcRecoLookupMapS::ProcessPhotonHit ( )

private

Definition at line 243 of file PndDrcRecoLookupMapS.cxx.

Referenced by Exec().

 {
   //Int_t nofChPho = 0;
   fDrcLutInfoArray->Clear();
   PndDrcLutInfo lutinfo;
   
   PndDrcPDPoint* Ppt=NULL;
   PndDrcPDHit* pdhit=NULL;
   PndMCTrack* tr = NULL;
   PndMCTrack* trMr;
         
   if (fVerbose > 0) {
     cout <<"PndDrcRecoLookupMapS: Number of Detected Photon MCPDPoints : "<<fPDPointArray->GetEntries()<<endl;
     cout <<"PndDrcRecoLookupMapS: Number of Detected Photon PDHits : "<<fPDHitArray->GetEntries()<<endl;
   }
   
  // Loop over PndDrcPDHits
  for(Int_t k=0; k<fPDHitArray->GetEntriesFast(); k++) {
       
     pdhit = (PndDrcPDHit*)fPDHitArray->At(k);
         
     Int_t mcPDRef= pdhit->GetRefIndex();
     
     if(mcPDRef < 0)continue;
     Ppt = (PndDrcPDPoint*)fPDPointArray->At(mcPDRef);
     
     if(Ppt->GetBarPointID() < 0) continue;    
     PndDrcBarPoint *fBarPoint= (PndDrcBarPoint*)fBarPointArray->At(Ppt->GetBarPointID());
     fBarId = fBarPoint->GetDetectorID();
     if(fBarId < 0)continue;
     //cout<<"bar name - "<<fGeoH->GetPath(fBarId)<<endl;
             
     Int_t trID= Ppt->GetTrackID();
     tr = (PndMCTrack*)fMCArray->At(trID);
     
     Int_t trMID=tr->GetMotherID();
     
     //if(trMID > -1){
     trMr = (PndMCTrack*)fMCArray->At(trMID);        
     Int_t trMpdg=trMr->GetPdgCode();
       
     if(trMr->GetMotherID()==-1){ // charged track is a primary
     
       lutinfo.SetChPartPdg(trMpdg);
     
       // expected cherenkov angle
       Double_t Mrmass;
       Int_t MoSign;
       if(fabs(trMpdg) == 11){Mrmass = 0.000511; MoSign = 1;}
       if(fabs(trMpdg) == 13){Mrmass = 0.105658;  MoSign = 1;}
       if(fabs(trMpdg) == 211){Mrmass = 0.139570; MoSign = -1;}
       if(fabs(trMpdg) == 321){Mrmass = 0.49368; MoSign = -1;}
       if(fabs(trMpdg) == 2212){Mrmass = 0.938; MoSign = -1;}
       if(fabs(trMpdg) == 50000050){Mrmass = 0.0; MoSign = -1;}
       Double_t Mrmom = trMr->GetMomentum().Mag();      
       CHexp = acos(sqrt(pow(Mrmom,2) + pow(Mrmass,2))/Mrmom/fGeo->nQuartz());
       lutinfo.SetCherenkovMC(CHexp);
       //cout<<"+++++++++++++++++++++++++++"<<endl;
       //cout<<"CH expected = "<<CHexp<<endl;
     
       fxPHit= pdhit->GetX();
       fyPHit= pdhit->GetY();
       fzPHit= pdhit->GetZ();
       ftime = pdhit->GetTime(); 
       //ftime = Ppt->GetTime();     
            
       //get the pixel id
       fpixID = pdhit->GetDetectorID();
       if(fpixID<0)continue;
        
       //cout<<"hit id "<<fpixID<<", x = "<<pdhit->GetX()<<", y = "<<pdhit->GetY()<<", time = "<<ftime<<endl;      
       
       Double_t xPPoi= Ppt->GetX();
       Double_t yPPoi= Ppt->GetY();
       
       lutinfo.AddTruePath(Ppt->GetLength());     
            
       // momentum of a photon on the PD Plane     
       fPphoPD.SetXYZ(Ppt->GetPx(), Ppt->GetPy(), Ppt->GetPz());
       
       {     
         // mother momentum/direction      
         fPMo.SetXYZ(trMr->GetMomentum().X(), trMr->GetMomentum().Y(), trMr->GetMomentum().Z()); 
         //fPMo.Print();
         lutinfo.SetChPartDir(fPMo);
         //cout<<"mother phi = "<<fPMo.Phi()/3.1415*180.<<endl;            
         if(fB > 0.){       
           Double_t PtMo = sqrt(pow(trMr->GetMomentum().X(),2)+pow(trMr->GetMomentum().Y(),2));
           Double_t Rratio = fR*fR/2./pow((PtMo/0.29979/fB)*100.,2);       
           Double_t phi_extra = fPMo.Phi() + MoSign * trMpdg/fabs(trMpdg) * acos(1. - Rratio);     
           fHAngleInBDeg = 0.5 * MoSign * trMpdg/fabs(trMpdg) * (acos(1. - Rratio) /TMath::Pi())*180.;
           //cout<<"B = "<<fB<<", R = "<<fR<<", trMpdg = "<<trMpdg<<", half angle in B = "<<fHAngleInBDeg<<endl;
           fPMo.SetPhi(phi_extra);
         }       
         lutinfo.SetChPartDirInBar2(fPMo);
       }
 
       {
         //mother momentum/direction from the bar point
         fBarPoint->Momentum(fPMo);
         lutinfo.SetChPartDirInBar(fPMo);        
       }
       
       //cout<<"Mother vector (global cs) : "<<endl;
       TVector3 motherMom = fPMo.Unit();
       //motherMom.Print();
          
       // initial momentum of the photon      
       fPphoInit = tr->GetMomentum();
       //cout<<"Initial momentum of the photon :"<<endl;
       //fPphoInit.Print();
                               
       // real = generated cherenkov angle
       CHreal = fPphoInit.Angle(fPMo);
       lutinfo.SetCherenkovReal(CHreal);
       //cout<<"CH real (generated) = "<<CHreal<<endl;
     
       Double_t etot = sqrt(pow(fPphoPD.X(),2) + pow(fPphoPD.Y(),2) +pow(fPphoPD.Z(),2));      
       flambdah=197.0*2.0*TMath::Pi()/(etot*1.0E9);//wavelength of photon in nm
       lutinfo.AddLambda(flambdah);
       
       // production point of the photon
       fStartVertex = tr->GetStartVertex();
       
       //cout<<"start: phi = "<<fStartVertex.Phi()<<endl;
       //fStartVertex.Print();
       ftime0 = tr->GetStartTime()/1.0E9;
       lutinfo.AddHitTime(ftime-ftime0);
       
       //cout<<"fstart time = "<<ftime0<<endl;
       //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
       // proceed only if photon was born in the bar for which the LUT is used!!!
       //if(fStartVertex.Phi() >0./180.*fpi && fStartVertex.Phi() < 6./180.*fpi){
       if(fPphoInit.Z() > 0. || fPphoInit.Z() < 0.){ // reflected photons
        
       // Transformation of the photon initial direction to the local bar coodr system:      
       // to get to the bar coord system.
       fPhiRot = InBarCoordSyst(fStartVertex, &fBBver1, &fBBver2, &fBBver3, &fBBver4);      
       
       // Photon initial momentum in the bar coord. syst.
       fPphoB = (fGeoH->MasterToLocalShortId(fPphoInit, fBarId) - fGeoH->MasterToLocalShortId((0.,0.,0.),fBarId)).Unit(); // vector
       //fPphoB.Print();      
       
       // mother momentum in the bar' (and bar) coord syst:
       fPMoB = (fGeoH->MasterToLocalShortId(fPMo, fBarId) - fGeoH->MasterToLocalShortId((0.,0.,0.),fBarId)).Unit(); // vector       
       //PMBar.Print();     
            
       // to find Cherenkov Phi (NOT THETA!!!):
       TVector3 kB;
       kB.SetXYZ(fkxBar, fkyBar, fkzBar);
             
       // use lookups to get the kXbar and kYbar of photons:
       NHitPix = fDigiPar->GetNHitPixels();
       //cout<<"N pixels = "<<NHitPix<<endl;
       NAmb = fDigiPar->GetNAmbiguities();
       //cout<<"N amb = "<<NAmb<<endl;
       NPixPar = fDigiPar->GetNPixelParam();
       Double_t par[NPixPar]; 
            
       if (fDigiPar->GetParamsForPixel(fpixID, par) == kFALSE) {
         cout<<"NO SUCH HIT IN LUT!!!"<<endl;
         continue;
       }
       
       //cout<<"LOOKUP!!!!"<<endl;      
       //cout<<"N Hit Pix = "<<NHitPix<<endl;
       //cout<<"N Amb = "<<NAmb<<endl;
       //cout<<"NPixPar = "<<NPixPar<<endl;
 /*      cout<<"par 0 = "<<par[0]<<endl; // kx 
       cout<<"par 1 = "<<par[1]<<endl; // ky 
       cout<<"par 2 = "<<par[2]<<endl; // kxB 
       cout<<"par 3 = "<<par[3]<<endl; // kyB 
       cout<<"par 4 = "<<par[4]<<endl; // kxU1
       cout<<"par 5 = "<<par[5]<<endl; // kyU1
       cout<<"par 6 = "<<par[6]<<endl; // kxU2
       cout<<"par 7 = "<<par[7]<<endl; // kyU2
       cout<<"par 8 = "<<par[8]<<endl; // kxU3
       cout<<"par 9 = "<<par[9]<<endl; // kyU3
       cout<<"par 10 = "<<par[10]<<endl; // kxBU1
       cout<<"par 11 = "<<par[11]<<endl; // kyBU1
       cout<<"par 12 = "<<par[12]<<endl; // kxBU2
       cout<<"par 13 = "<<par[13]<<endl; // kyBU2
       cout<<"par 14 = "<<par[14]<<endl; // kxBU3
       cout<<"par 15 = "<<par[15]<<endl; // kyBU3
       cout<<"par 16 = "<<par[16]<<endl; // kxUB 
       cout<<"par 17 = "<<par[17]<<endl; // kyUB 
       cout<<"par 18= "<<par[18]<<endl; // kxUU 
       cout<<"par 19= "<<par[19]<<endl; // kyUU 
       cout<<"par 20= "<<par[20]<<endl; // kxUUU
       cout<<"par 21= "<<par[21]<<endl; // kyUUU
       cout<<"par 22= "<<par[22]<<endl; // kxBUU
       cout<<"par 23= "<<par[23]<<endl; // kyBUU
       cout<<"par 24= "<<par[24]<<endl; // kxBUB 
       cout<<"par 25= "<<par[25]<<endl; // kyBUB 
       cout<<"par 26= "<<par[26]<<endl; // kxUBU
       cout<<"par 27= "<<par[27]<<endl; // kyUBU         
 */     
       // fill in NPixPar/2 * 1 ambiguities:
       Double_t kX, kY, kZ;                                  
       std::vector<TVector3> ambig;
       std::vector<Double_t> CHreco;     
       std::vector<Double_t> Tamb; 
       std::vector<Double_t> Path; 
                 
       // length of the vector should be the same as NumAmb*8   !!!!!
       Int_t NumAmb = NAmb;  //=14        
       TVectorD CHdiff(8*NumAmb); 
       TVectorD Adiff(8*NumAmb); // vector to store angle diff btw real and reco photons
             
       for(Int_t i=0; i<NumAmb; i++){
         // in bar' coordinate system    
         kX = par[i*2]; // kXD   
         kY = par[i*2 + 1]; // kYD       
         kZ = -sqrt(1. - pow(kX,2.) - pow(kY,2.));       
         //cout<<"mom:  x = "<<fPMoB.X()<<",  y = "<<fPMoB.Y()<<", z = "<<fPMoB.Z()<<endl;       
         fkBar.SetXYZ(kX, kY, kZ);
         //cout<<"photon from LUT in bar' coordinate system: "<<endl;
         //fkBar.Print();
                 
         for(Int_t jamb=0; jamb <8; jamb++){
                 if(jamb == 1){fkBar.SetXYZ(-kX, kY, kZ);}
                 if(jamb == 2){fkBar.SetXYZ(kX, -kY, kZ);}
                 if(jamb == 3){fkBar.SetXYZ(kX, kY, -kZ);}
                 if(jamb == 4){fkBar.SetXYZ(-kX, -kY, kZ);}
                 if(jamb == 5){fkBar.SetXYZ(kX, -kY, -kZ);}
                 if(jamb == 6){fkBar.SetXYZ(-kX, kY, -kZ);}
                 if(jamb == 7){fkBar.SetXYZ(-kX, -kY, -kZ);}
                                                         
                 // check if there is such enrty in the LUT and if the ambiguity fulfills the total internal reflection requitrement
                 if((kX != 0. && kY != 0.) && ((fkBar.Cross(fnX1)).Mag() > 1.00028/fGeo->nQuartz() || 
                    (fkBar.Cross(fnY1)).Mag() > 1.00028/fGeo->nQuartz())){
                         ambig.push_back(fkBar);
                         CHreco.push_back(fkBar.Angle(fPMoB));
                         CHdiff(8*i+jamb) = fabs(CHreco[8*i+jamb] - CHexp);
                         Tamb.push_back(RecoAmbigTime(fkBar, fStartVertex, &fPath, 0));                  
                         Path.push_back(fPath);   
                         Adiff(8*i+jamb) = fkBar.Angle(fPphoB);
                         //cout<<"reco ch angle = "<<fkBar.Angle(fPMoB)<<endl;
                         //cout<<"number of bounces = "<<NumberOfBounces(fStartVertex, fPphoB, fBarId)<<endl;      
                         //cout<<"CHdiff = "<<CHdiff(8*i+jamb)<<", CHreco"<<8*i+jamb<<" = "<<CHreco[8*i+jamb]<<endl;
                         
                         //fill   only with credible information:
                         lutinfo.AddAngle(fkBar.Angle(fPMoB));
                         lutinfo.AddTime(RecoAmbigTime(fkBar, fStartVertex, &fPath, 0));
                         lutinfo.AddPath(fPath);
                         lutinfo.AddChDiff(CHreco[8*i+jamb] - CHexp);
                         //cout<<"fkBar = "<<endl;
                         //fkBar.Print();
                         lutinfo.AddNOfBounces(NumberOfBounces(fStartVertex, fkBar/*fPphoB*/, fBarId));
                         
                         //fkBarXHist->Fill(kX, fPphoB.X()); 
                         //fkBarYHist->Fill(kY, fPphoB.Y());                               
                 }else{
                         ambig.push_back((0.,0.,0.));
                         CHreco.push_back(-999.);
                         CHdiff(8*i+jamb) = -999.;
                         Tamb.push_back(-999.);                  
                         Path.push_back(-999.);   
                         Adiff(8*i+jamb) = -999.;                                        
                 }
         }// for                         
       }
                 
         fPhiMap = TMath::Nint(fPMo.Phi()/fpi*180.);
         fThetaMap = TMath::Nint(fPMo.Theta()/fpi*180.);
         //cout<<"phi = "<<fPMo.Phi()/fpi*180.<<", phi map = "<<fPhiMap<<", theta map = "<<fThetaMap<<endl;
           
         if(fB > 0.){
          // fMapHist->Fill(fThetaMap,fPhiMap-fHAngleInBDeg);  // assuming that R = 50 cm  
         }            
         if(fB == 0.){
          // fMapHist->Fill(fThetaMap,fPhiMap);    
         }                   
     } // if start vertex lies within the right bar!
     }// photon from primary particle 
              
  }// photon hits 
  new ((*fDrcLutInfoArray)[fDrcLutInfoArray->GetEntriesFast()]) PndDrcLutInfo(lutinfo);
 }

void PndDrcRecoLookupMapS::ProcessPhotonMC ( )

private

Definition at line 213 of file PndDrcRecoLookupMapS.cxx.

References fMCArray, and fVerbose.

 {
   if (fVerbose > 0) {
     cout <<"PndDrcRecoLookupMapS: Number of Detected MC Tracks : "<<fMCArray->GetEntries()<<endl;
     PndMCTrack* ptr = NULL;
     PndMCTrack* trMr=NULL;
     
   }
 }

Double_t PndDrcRecoLookupMapS::RecoAmbigTime	(	TVector3	kb,
		TVector3	start,
		Double_t *	l,
		Bool_t	printout
	)

private

Definition at line 644 of file PndDrcRecoLookupMapS.cxx.

References PndGeoDrc::barBoxZDown(), PndGeoDrc::barBoxZUp(), cos(), Double_t, PndGeoDrc::EVlen(), fabs(), fGeo, fpi, L1, PndGeoDrc::nEV(), PndGeoDrc::nQuartz(), and sin().

Referenced by ProcessPhotonHit().

                                                                                                      {  // [ns]  
 
   // group velocities for oil and quartz for Noil and Nquartz assuming that mean lambda = 410 nm:
   //Double_t nOil = fGeo->nEV();
   Double_t u_oil = 19.8;//30./nOil;
   //Double_t u_quartz = 30./1.46907; // for lambda = 410 nm 
   //Double_t u_quartz = 30./1.47012; // for lambda = 400 nm 
   //Double_t u_quartz = 30./1.46838; // for lambda = 417 nm 
   //Double_t u_quartz = 30./1.47125; // for lambda = 390 nm 
   //Double_t u_quartz = 30./1.47248; // for lambda = 380 nm
   Double_t u_quartz = 20.;//30./fNquartz;//1.47805; // for lambda = 400 nm
   
   // kb - photon momentum at the production point, first calculate path in QUARTZ BAR:
   // assume that there are only DIRECT photons!!!   
   
   // path in quartz:
   Double_t Z0 = 0.;
   Double_t kz = 0.;
   if(kb.Z() <  0.){Z0 = start.Z() - fGeo->barBoxZUp(); kz = kb.Z();}
   if(kb.Z() >= 0.){Z0 = 2.*fGeo->barBoxZDown() - fGeo->barBoxZUp() - start.Z(); kz = -kb.Z();}
   
   if(printout){
     cout<<"-I- Reco ambig Time: Zstart = "<<start.Z()<<", Z0 = "<<Z0<<endl;
   }
   
   Double_t L0 = Z0/fabs(cos(kb.Theta()));
   kb.SetZ(kz);
   Double_t thetaOil = asin(fGeo->nQuartz()/fGeo->nEV()*sin(kb.Theta()));
   if(printout){
     cout<<"RECO ambig TIME: thetas are: quartz = "<<kb.Theta()/fpi*180.<<", oil = "<<thetaOil/fpi*180.<<endl;
   }
   kb.SetTheta(thetaOil);
   Double_t Z1 = fGeo->EVlen();
   Double_t L1 = Z1/fabs(cos(kb.Theta())); // path in oil
     
   if(printout){
     cout<<"-I- Reco ambig Time: Zev = "<<fGeo->EVlen()<<endl;     
     cout<<"-I- Reco ambig Time: L0 = "<<L0<<", L1 = "<<L1<<endl;
     cout<<"-I- Reco ambig Time: time is "<< L0/u_quartz + L1/u_oil <<endl; // 19.5 cm/ns - average (lambda) group speed
   } 
   
   *l = L0 + L1;
    
   return L0/u_quartz + L1/u_oil; //[ns]  
 }

void PndDrcRecoLookupMapS::SetDefaultParameters ( )

private

Set the parameters to the default values.

void PndDrcRecoLookupMapS::SetHistoWidth ( double wid = 0.01 )

inline

Definition at line 70 of file PndDrcRecoLookupMapS.h.

References fWidth.

70 {fWidth = wid;} // 0.01

PndDrcRecoLookupMapS::fWidth

Double_t fWidth

Definition: PndDrcRecoLookupMapS.h:141

void PndDrcRecoLookupMapS::SetMagneticField ( double mag = 0. )

inline

Definition at line 72 of file PndDrcRecoLookupMapS.h.

References fB.

72 {fB = mag;} // [Tesla]

PndDrcRecoLookupMapS::fB

Double_t fB

Definition: PndDrcRecoLookupMapS.h:148

void PndDrcRecoLookupMapS::SetNSigma ( double nsig = 3. )

inline

Definition at line 71 of file PndDrcRecoLookupMapS.h.

References fNSigma, and nsig.

71 {fNSigma = nsig;}

PndDrcRecoLookupMapS::fNSigma

Double_t fNSigma

Definition: PndDrcRecoLookupMapS.h:147

nsig

int nsig

Definition: toy_core.C:46

void PndDrcRecoLookupMapS::SetOutputFile ( TString outName = "PndDrcRecoLookupMapS_output.root" )

inline

Definition at line 73 of file PndDrcRecoLookupMapS.h.

References fOutputName.

73 {fOutputName = outName;}

PndDrcRecoLookupMapS::fOutputName

TString fOutputName

Definition: PndDrcRecoLookupMapS.h:140

void PndDrcRecoLookupMapS::SetParContainers ( )

privatevirtual

Definition at line 172 of file PndDrcRecoLookupMapS.cxx.

References fDigiPar, fGeoH, fVerbose, PndDrcDigiPar::GetNAmbiguities(), PndGeoHandling::Instance(), run, and PndGeoHandling::SetParContainers().

                                             {
 
   // Get run and runtime database
   FairRunAna* run = FairRunAna::Instance();
   if ( ! run ) Fatal("SetParContainers", "No analysis run");
 
   FairRuntimeDb* db = run->GetRuntimeDb();
   if ( ! db ) Fatal("SetParContainers", "No runtime database");
 
   // Get DIRC digitisation parameter container
   fDigiPar = (PndDrcDigiPar*)(db->getContainer("DIRCLookupTable"));
   cout << "-I- PndDrcRecoLookupMapS::SetParContainers(). read parameters" << endl;
   cout<<"read them!"<<endl;
   //cout << "-I- PndDrcRecoLookupMapS: Number of hit pixels "<< fDigiPar->GetNHitPixels()<< endl;
   cout << "-I- PndDrcRecoLookupMapS: Number of hit pixels "<< fDigiPar->GetNAmbiguities()<< endl;
   
   if ( fGeoH == NULL ) fGeoH = PndGeoHandling::Instance();
   fGeoH->SetParContainers(); 
   if(fVerbose>1) Info("SetParContainers","done.");
   return;
 }