#include <PndRichRecoTask.h>

Inheritance diagram for PndRichRecoTask:

Public Member Functions
	PndRichRecoTask ()

	~PndRichRecoTask ()

virtual InitStatus	Init ()

virtual void	Exec (Option_t *opt)

void	SetPersistence (Bool_t persistence)

void	SetGeoVersion (UInt_t version)

void	SetParticleID (UInt_t pid)

void	SetTrackPosition (TVector3 pos)

void	SetTrackDirection (TVector3 dir)

void	SetTrackPositionSecond (TVector3 pos)

void	SetTrackDirectionSecond (TVector3 dir)

void	SetNumberOfEvents (UInt_t events)

TF1 *	GetPeakParameters (std::vector< Double_t > v, UInt_t nbins, Double_t xmin, Double_t xmax, Double_t sig)

void	FinishEvent ()

void	FinishTask ()

Private Member Functions
	ClassDef (PndRichRecoTask, 1)
	PndRichBarPoint TCA. More...

Private Attributes
Bool_t	fPersistence

UInt_t	fNumberOfEvents

UInt_t	fGeoVersion

UInt_t	fParticleID

UInt_t	fEvent

PndRichReco *	fRichReco

TVector3	fTrackPosition

TVector3	fTrackDirection

TVector3	fTrackPositionSecond

TVector3	fTrackDirectionSecond

std::vector< UInt_t >	vnhits

std::vector< Double_t >	vth

std::vector< Double_t >	vph

std::vector< Double_t >	vdth

std::vector< Double_t >	vmean

std::vector< Double_t >	vsigma

PndRichResolution *	fRichResolution

TClonesArray *	fRichBarPoint

Detailed Description

Definition at line 15 of file PndRichRecoTask.h.

Constructor & Destructor Documentation

PndRichRecoTask::PndRichRecoTask ( )

Default constructor

Definition at line 22 of file PndRichRecoTask.cxx.

References fVerbose.

                                  :
   FairTask("Rich Reco task")  ,
   fPersistence(kTRUE),
   fNumberOfEvents(0),
   fEvent(0),
   fTrackPositionSecond(0,0,0),
   fTrackDirectionSecond(0,0,0)
 {  
   fVerbose=1;
 }

PndRichRecoTask::~PndRichRecoTask ( )

Destructor

Definition at line 35 of file PndRichRecoTask.cxx.

35 {

36 }

Member Function Documentation

PndRichRecoTask::ClassDef	(	PndRichRecoTask	,
		1
	)

private

PndRichBarPoint TCA.

void PndRichRecoTask::Exec ( Option_t * opt )

virtual

Virtual method Exec

Definition at line 66 of file PndRichRecoTask.cxx.

References dbpoint::beta, Double_t, fEvent, fRichBarPoint, fRichReco, fRichResolution, fTrackDirection, fTrackDirectionSecond, fTrackPosition, fTrackPositionSecond, fVerbose, PndRichReco::GetDThetas(), PndRichBarPoint::GetPdgCode(), PndRichReco::GetPhis(), PndRichReco::GetThetas(), i, dbpoint::mass, mom, p, chigen::pythia::pdf, dbpoint::phi, pnt, pos, PndRichReco::RichFullReconstruction(), PndRichResolution::Shift(), PndRichResolution::Sigma(), CAMath::Sqrt(), dbpoint::theta, vdth, vmean, vnhits, vph, vsigma, vth, dbpoint::x, and dbpoint::y.

                                     {
   if(fVerbose > 0) {
     //cout << "==================== EVENT " << evt << endl;
   }
    fEvent++;
    //fRichReco->RichFullReconstruction();
    
    TVector3 pos = fTrackPosition;
    TVector3 dir = fTrackDirection.Unit();
    Double_t mom = 0;
    Double_t tm = 0;
    Int_t pdg = 0;
    
    PndRichBarPoint *richHit = NULL;
    Int_t richEntries = fRichBarPoint->GetEntriesFast();
    if (richEntries) {
       richHit = (PndRichBarPoint*)fRichBarPoint->At(0);
       //pos = richHit->GetPosition0();
       //dir = richHit->GetMomentum0().Unit();
       //mom = richHit->GetMomentum0().Mag();
       pos = TVector3(richHit->GetX(),richHit->GetY(),richHit->GetZ());
       TVector3 p = TVector3(richHit->GetPx(),richHit->GetPy(),richHit->GetPz());
       dir = p.Unit();
       mom = p.Mag();
       tm = richHit->GetTime();
       pdg = richHit->GetPdgCode();
    }
    
    Int_t richPhot = 0;
    Float_t richThetaC = -1000, richThetaCErr = 0;
    Float_t richQuality = 1000000;
    // first particle
    fRichReco->RichFullReconstruction(pos,dir,tm,
                                      richQuality,richThetaC,richThetaCErr,richPhot);
    std::vector<Double_t> dth = fRichReco->GetDThetas();
 
    std::cout << std::setprecision(10) << std::endl;
    std::cout << "PndRichRecoTask: " << fEvent << " " << richQuality << " " << richThetaC
       << " " << richThetaCErr << " " << richPhot << " " << 1
       << " " << mom << " " << dir.Phi() << " " << dir.Theta()
       << " " << pos.X() << " " << pos.Y() << " " << pdg << std::endl;
    Double_t mass[5] = {0.000510998961,0.1056583745,0.13957018,0.493667,0.9382720813};
    Double_t pdf[5];
    Double_t pdfsum = 0;
    dbpoint pnt;
    for(Int_t i=0;i<5;i++)
    {
       Float_t richThetaCx = richThetaC;
       pnt.mass = mass[i];
       pnt.beta = richThetaCx;
       pnt.x = pos.X();
       pnt.y = pos.Y();
       pnt.theta = dir.Theta();
       pnt.phi = dir.Phi();
       pnt.beta += fRichResolution->Shift(pnt);
       richThetaCx += fRichResolution->Shift(pnt);
       
       std::cout << i << " " << fRichResolution->Shift(pnt) << " " << fRichResolution->Sigma(pnt) << " " << mass[i] << " " <<
          pnt.beta << " " << pnt.x << " " << pnt.y << " " << pnt.theta << " " << pnt.phi << std::endl;
       
       Double_t beta = mom / TMath::Sqrt(mom*mom + mass[i]*mass[i]);
       TF1 *gausPdf = new TF1("gausPdf","gausn",0,1);
       gausPdf->SetParameter(0,1);
       gausPdf->SetParameter(1,richThetaCx);
       gausPdf->SetParameter(2,richThetaCErr);
       //gausPdf->SetParameter(2,fRichResolution->Sigma(pnt));
       pdf[i] = gausPdf->Eval(beta);
       pdfsum += pdf[i];
       delete gausPdf;
    }
    std::cout << "particles PID: ";
    for(Int_t i=0;i<5;i++)
       std::cout << pdf[i] << " " << (pdfsum?pdf[i]/pdfsum:0.2) << " " ;
    std::cout << std::endl;
    
    vdth.insert(vdth.end(),dth.begin(),dth.end());
    vnhits.push_back(dth.size());
    vmean.push_back(richThetaC);
    vsigma.push_back(richThetaCErr);
    
    std::vector<Double_t> th = fRichReco->GetThetas();
    std::vector<Double_t> ph = fRichReco->GetPhis();
    vth.insert(vth.end(),th.begin(),th.end());
    vph.insert(vph.end(),ph.begin(),ph.end());
 
    for(size_t i=0; i<ph.size(); i++)
       std::cout << "::HitPars()  " << fEvent << " " << ph.at(i) << " " <<
          th.at(i) << " " <<  dth.at(i) << std::endl;
    
    if (fTrackDirectionSecond.Mag()!=0) {
       // secon particle
       fRichReco->RichFullReconstruction(fTrackPositionSecond,fTrackDirectionSecond.Unit(),0.,
                                         richQuality,richThetaC,richThetaCErr,richPhot);
       std::cout << "PndRichRecoTask: " << fEvent << " " << richQuality << " " << richThetaC
          << " " << richThetaCErr << " " << richPhot << " " << 2 << std::endl;
    }
 }

void PndRichRecoTask::FinishEvent ( )

Definition at line 164 of file PndRichRecoTask.cxx.

165 {

166 }

void PndRichRecoTask::FinishTask ( )

Definition at line 168 of file PndRichRecoTask.cxx.

References Double_t, f, f1, fabs(), GetPeakParameters(), i, mean, nhits, vdth, vmean, and vsigma.

 {
    std::cout << "::FinishTask()  " << vmean.size() << " " << vsigma.size() << std::endl;
    TF1 * f = GetPeakParameters(vmean,1100,0.95,1.005,0.0003);
    
    double chisq=f->GetChisquare();
    double ndf=f->GetNDF();
    double chisqdf=chisq/ndf;
    std::cout << "Chisquare: " << chisq << "/" << ndf << " : " << chisqdf << std::endl;
    
    //Double_t amp = f->GetParameter(0); //value of 0th parameter //[R.K. 03/2017] unused variable?
    //Double_t eamp = f->GetParError(0); //error on 0th parameter //[R.K. 03/2017] unused variable?
    Double_t mean = f->GetParameter(1); //value of 1st parameter
    Double_t emean = f->GetParError(1); //error on 1st parameter
    Double_t sig = f->GetParameter(2); //value of 1st parameter
    Double_t esig = f->GetParError(2); //error on 1st parameter
 
    size_t nhits = 0;
    for(size_t i=0; i<vmean.size(); i++) {
       if ( std::fabs(vmean.at(i)-mean) < 3*sig ) nhits++;
    }
 
    std::cout << "::FinishTask() " << mean << " " << sig << std::endl;
    std::cout << "CalData: " << chisq << std::endl;
    std::cout << "CalData: " << ndf << std::endl;   
    std::cout << "CalData: " << nhits << std::endl;
    std::cout << "CalData: " << mean << std::endl;
    std::cout << "CalData: " << emean << std::endl;
    std::cout << "CalData: " << sig << std::endl;
    std::cout << "CalData: " << esig << std::endl;
    
    TF1 * f1 = GetPeakParameters(vdth,500,-0.05,0.05,0.003);
    
    chisq=f1->GetChisquare();
    ndf=f1->GetNDF();
    chisqdf=chisq/ndf;
    std::cout << "Chisquare: " << chisq << "/" << ndf << " : " << chisqdf << std::endl;
    
    //amp = f1->GetParameter(0); //value of 0th parameter //[R.K. 03/2017] unused variable?
    //eamp = f1->GetParError(0); //error on 0th parameter //[R.K. 03/2017] unused variable?
    mean = f1->GetParameter(1); //value of 1st parameter
    emean = f1->GetParError(1); //error on 1st parameter
    sig = f1->GetParameter(2); //value of 1st parameter
    esig = f1->GetParError(2); //error on 1st parameter
 
    std::cout << "::FinishTask() " << mean << " " << sig << std::endl;
    std::cout << "CalData: " << chisq << std::endl;
    std::cout << "CalData: " << ndf << std::endl;   
    std::cout << "CalData: " << mean << std::endl;
    std::cout << "CalData: " << emean << std::endl;
    std::cout << "CalData: " << sig << std::endl;
    std::cout << "CalData: " << esig << std::endl;
    
 }

TF1 * PndRichRecoTask::GetPeakParameters	(	std::vector< Double_t >	v,
		UInt_t	nbins,
		Double_t	xmin,
		Double_t	xmax,
		Double_t	sig
	)

Definition at line 223 of file PndRichRecoTask.cxx.

References Double_t, i, and xmin.

Referenced by FinishTask().

 {
    TH1F * hbetafit = new TH1F("hist", "for gauss fit", nbins, xmin, xmax);
    for(size_t i=0; i<v.size(); i++)
       hbetafit->Fill(v.at(i));
    int ipeak = hbetafit->GetMaximumBin();
    double amp0 = hbetafit->GetMaximum();
    double peak = xmin + (ipeak-0.5)*(xmax-xmin)/nbins;
    TF1 * gfit = new TF1("Gaussian","gaus",xmin,xmax);
    gfit->SetParameter(0,amp0);
    gfit->SetParameter(1,peak);
    gfit->SetParameter(2,sig);
    for(size_t i=0; i<5; i++) {
       Double_t meanf = gfit->GetParameter(1);
       Double_t sigf = gfit->GetParameter(2);
       double xminf = meanf - 3*sigf;
       double xmaxf = meanf + 3*sigf;
       hbetafit->Fit("Gaussian","Q","",xminf,xmaxf);
    }
    delete hbetafit;
    return gfit;
 }

InitStatus PndRichRecoTask::Init ( )

virtual

Virtual method Init

Definition at line 40 of file PndRichRecoTask.cxx.

References fGeoVersion, fNumberOfEvents, fRichBarPoint, fRichReco, fRichResolution, vmean, vnhits, and vsigma.

                                  {
 
   FairRootManager *ioman = FairRootManager::Instance();
   
   if (!ioman) 
    {
       cout << "-E- PndRichRecoTask: "
          << "RootManager not instantised!" << endl;
       return kFATAL;
    }
 
    fRichReco = new PndRichReco(fGeoVersion); // 13 - one of the flat mirror variant
    vnhits.reserve(fNumberOfEvents);
    vmean.reserve(fNumberOfEvents);
    vsigma.reserve(fNumberOfEvents);
    fRichBarPoint = dynamic_cast<TClonesArray *> (ioman->GetObject("RichBarPoint"));
 
    fRichResolution = new PndRichResolution();
 
   cout << "-I- PndRichRecoTask: Intialisation successfull " << endl;
   return kSUCCESS;
 }