// $Id: TrackDirection2Module.cc,v 1.16 2005/10/21 13:30:08 moreggia Exp $
// Author: elena Nov, 19 2004
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: TrackDirection2Module *
* Package: Fitting *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// TrackDirection2Module
//
// This module is devoted to the reconstruction of the shower direction.
// It implements some different algorithms (descibed in detail in ....)
// Essentially the module keeps the points found by pattern recogniton
// (clustering or Hough transform) and first find the plane that contains the
// track and the detector (TDP).
// The reserarch of TDP can be done in the following ways:
// - further selection of points with hough transform
// - further selection of points with shape selection method
// - no further selection
// In each case the TDP is founded by a fit (least squares, median or hough)
// of the x-t, y-t projections of points on the plane z=0.
//
// Then the shower direction is reconstructed by one of the following methods:
// - analytical approximated 1 AA1()
// - analytical approximated 2 AA2()
// - numerical exact 1 NE1()
// - numerical exact 2 NE2()
// - analytical exact 1 AE1()
//
// The AA1() method is in each case used in order to initialize the fit for all
// other methods.
//
// Config file parameters
// ======================
//
// fNumHitsMinimum : minimum number of hits per pixel in a GTU
// fNumPointsMinimum : minimum number of points to proceed with reconstruction
// fUseHough [bool] : use hough trasnform to find the TDP
// fUseHoughwithselection [bool] : use the points selected with Hough transform
// in the reconstruction of direction
// fUseShape [bool] : use shape selection method to find the TDP
// fUseShapeInModules [bool] : use the points selected with shape selction
// in the reconstruction of direction
// fDebugInfo [bool] : compute and display some debug informations
//
// fMulti1 : multiplier for shape selection method
// fMulti2 : multiplier for shape selection method
//
// fAA1FitMethod : fit method for AA1 algorithm
// - Valid options : linear (least squares fit)
// median (median fit)
// hough (hough fit)
//
// fMethod : direction reconstruction method
// - Valid options : AA1 || AA2 || NE1 || NE2 || AE1 (single algorithm)
// all (executes all algorithms)
//
// fFixTmaxNumeric [bool] : fix shower maximum parameters in numerical fits
// fErrAngle : angular error [deg]
//
// fStat1 : angular error [deg] value 1 for statistics
// fStat2 : angular error [deg] value 2 for statistics
// - Statistics are made between 0 < err < fStat1 and fStat1 < err < fStat2
//
// fDoGraphUseShape [bool] : save some debug graph in rootfile
//
// fMinuitOutputLevel : set the MINUIT output display level
//
#include "EusoCluster.hh"
#include "TrackDirection2Module.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "LeastSquaresFit.hh"
#include "MedianFit.hh"
#include "HoughFit.hh"
#include "RecoRootEvent.hh"
#include "EConst.hh"
#include "Config.hh"
#include "ERunParameters.hh"
#include <TGraph.h>
#include <TCanvas.h>
#include <TMinuit.h>
#include <TVector3.h>
#include <TH2F.h>
#include <TLegend.h>
using namespace TMath;
using namespace sou;
using namespace EConst;
Double_t DEG = RadToDeg();
Double_t PI = Pi();
Double_t RT = EarthRadius()/km; //earth radius in km
Double_t CLUCE = Clight()*1.e-3; //constant c in unit km/microsecond
//FCN function called by numeric method NE1
void ChiSquareNE1(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag);
//FCN function called by numeric method NE2
void ChiSquareNE2(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag);
ClassImp(TrackDirection2Module)
ClassImp(ContainerData)
//_____________________________________________________________________________
TrackDirection2Module::TrackDirection2Module() : RecoModule("TrackDirection2") {
//
// ctor
//
}
//_____________________________________________________________________________
TrackDirection2Module::~TrackDirection2Module() {
//
// dtor
//
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::Init() {
//
// Initialization of variables
//
Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
useLFaa1 = kFALSE; useMFaa1 = kFALSE; useHFaa1 = kFALSE;
useLFplane = kFALSE; useMFplane = kFALSE; useHFplane = kFALSE;
fStat1 = Conf()->GetNum("TrackDirection2Module.fStat1");
fStat2 = Conf()->GetNum("TrackDirection2Module.fStat2");
fDoGraphUseShape = Conf()->GetBool("TrackDirection2Module.fDoGraphUseShape");
fNumPointsMin = (Int_t)Conf()->GetNum("TrackDirection2Module.fNumPointsMin");
fNumHitsMinimum = (Int_t)Conf()->GetNum("TrackDirection2Module.fNumHitsMinimum");
fErrAngle = Conf()->GetNum("TrackDirection2Module.fErrAngle")/DEG;
if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="linear") useLFplane = kTRUE;
else if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="median") useMFplane = kTRUE;
else if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="hough") useHFplane = kTRUE;
else Msg(EsafMsg::Panic) << "Wrong config fTDPFitMethod value in TrackDirection2Module.cfg" << MsgDispatch;
if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="linear") useLFaa1 = kTRUE;
else if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="median") useMFaa1 = kTRUE;
else if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="hough") useHFaa1 = kTRUE;
else Msg(EsafMsg::Panic) << "Wrong config fAA1FitMethod value in TrackDirection2Module.cfg" << MsgDispatch;
if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AA1") fMethodIdentifier = 1;
else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="NE1") fMethodIdentifier = 2;
else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AA2") fMethodIdentifier = 3;
else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="NE2") fMethodIdentifier = 4;
else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AE1") fMethodIdentifier = 5;
else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="all") fMethodIdentifier = 6;
else Msg(EsafMsg::Panic) << "Wrong config fMethod value in TrackDirection2Module.cfg" << MsgDispatch;
fDoHough = Conf()->GetBool("TrackDirection2Module.fUseHough");
fOptionSelectionHough = Conf()->GetBool("TrackDirection2Module.fUseHoughwithselection");
fDoShapeSelection = Conf()->GetBool("TrackDirection2Module.fUseShape");
if ( fDoShapeSelection) {
fMulti1 = Conf()->GetNum("TrackDirection2Module.fMulti1");
fMulti2 = Conf()->GetNum("TrackDirection2Module.fMulti2");
}
fUseShapeSelectioninModules = Conf()->GetBool("TrackDirection2Module.fUseShapeInModules");
fFixTmaxNumeric = Conf()->GetBool("TrackDirection2Module.fFixTmaxNumeric");
fDebugInfo = Conf()->GetBool("TrackDirection2Module.fDebugInfo");
fMinuitOutputLevel = (Int_t)Conf()->GetNum("TrackDirection2Module.fMinuitOutputLevel");
if ( fMinuitOutputLevel < -1 || fMinuitOutputLevel > 3)
Msg(EsafMsg::Panic) << "Wrong fMinuitOutputLevel value in TrackDirection2Module.cfg" << MsgDispatch;
ConfigFileParser *pConfigGeneralEuso = Config::Get()->GetCF("General","Euso");
fHISS = pConfigGeneralEuso->GetNum("Euso.fAltitude");
// gtu length in microseconds
ConfigFileParser *pConfigMacroCell = Config::Get()->GetCF("Electronics","MacroCell");
fGtuLength = pConfigMacroCell->GetNum("MacroCell.fGtuTimeLength")/microsecond;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::PreProcess() {
//
// Pre-process of the reco event
//
fData.Clear();
fAA1done = kFALSE;
fAA1nan = kFALSE;
fNumPoints = 0;
fNumHits = 0;
fNumPointsSel = 0;
fNumHitsSel = 0;
fQuality = -1;
fCentroid.SetXYZ(0,0,0);
fNorm.SetXYZ(0,0,0);
fNormsel.SetXYZ(0,0,0);
fW.SetXYZ(0,0,0);
fU.SetXYZ(0,0,0);
fTrueDir.SetXYZ(0,0,0);
fTrueNorm.SetXYZ(0,0,0);
fTrueMax.SetXYZ(0,0,0);
fEASDir.SetXYZ(0,0,0);
fAngularSpeed = 0;
fBeta = 0;
fBetaInit = 0;
fHmax = 0;
fRmax = 0;
fTrueTheta = 0;
fTruePhi = 0;
fTHETAloc = 0;
fTHETAreco = 0;
fPHIreco = 0;
fTmaxFit = 0;
fDeltaTheta = 0;
fDeltaPhi = 0;
fDeltaEASDir = 0;
fDeltaTDP = 0;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::Process(RecoEvent *ev) {
fEv = ev;
RecoEventHeader header = fEv->GetHeader();
ERunParameters *runpars = header.GetRunPars();
if(fDoGraphUseShape){
gDirectory->cd("/");
string pathname = Form("TD2ProcessRecoEvent%d", fEv->GetHeader().GetNum());
TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
path->cd();
}
const RecoModuleData *gcm = fEv->GetModuleData("GTUClustering");
const RecoModuleData *moduledataHough = fEv->GetModuleData("HoughTransform");
if ( gcm == NULL && moduledataHough == NULL ) {
Msg(EsafMsg::Warning) << "Not found GTUClusteringModule and not found HoughTransformModule" << MsgDispatch;
Int_t fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
Msg(EsafMsg::Info) <<"fTotalNumPoints = "<<fTotalNumPoints<<MsgDispatch;
if ( fTotalNumPoints <= 0 ) {
Msg(EsafMsg::Warning) << "No points in event! Exit from this event." << MsgDispatch;
return kFALSE;
}
fPointsId.clear();
fNumPoints = 0;
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
if ( fEv->GetRecoPixelData(i)->GetCounts() >= fNumHitsMinimum ) {
fPointsId.push_back(i);
fNumPoints++;
}
}
Msg(EsafMsg::Info) << "Num. of pixels with at least " << fNumHitsMinimum << " hits = " << fNumPoints << MsgDispatch;
} else if(! (gcm == NULL) ) {
if ( gcm->GetObj("CluPixels") == NULL ) {
Msg(EsafMsg::Warning) << "No data in event" << MsgDispatch;
return kFALSE;
}
fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
if (fPointsId.size() == 0) {
fNumPoints = 0;
Msg(EsafMsg::Warning) << "No pixels selected: terminated." << MsgDispatch;
return kFALSE;
} else {
Msg(EsafMsg::Info) << "Pattern Recognition GtuClusteringModule data found with " << fPointsId.size()<< " points" << MsgDispatch;
fNumPoints = fPointsId.size();
}
} else if ( !(moduledataHough->GetObj("Clusters") == NULL) ) {
vector<EusoCluster*> cluHough;
cluHough= *(vector<EusoCluster*>*) moduledataHough->GetObj("Clusters");
fPointsId = cluHough[0]->GetPixelVector();
if (fPointsId.size() == 0) {
fNumPoints = 0;
Msg(EsafMsg::Warning) << "No data in event" << MsgDispatch;
return kFALSE;
} else {
Msg(EsafMsg::Info) << "Pattern Recognition HoughModule data found with " << fPointsId.size() << " points" << MsgDispatch;
fNumPoints = fPointsId.size();
}
}
fNumHits = 0;
for(Int_t i=0; i<fNumPoints; i++) {
RecoPixelData *pix = fEv->GetRecoPixelData( fPointsId[i] );
TVector3 pos(0,0,0);
Double_t phitmp = pix->GetPhi()+PI;
pos.SetMagThetaPhi( 1,pix->GetTheta(),phitmp );
fData.fSpPos.push_back( pos );
Int_t pixid = pix->GetPixelId();
fData.fPixXYZ.push_back( runpars->PixelCenter(pixid) );
fData.fSigmaPhi.push_back( pix->GetPhiSigma() );
fData.fSigmaTheta.push_back( pix->GetThetaSigma() );
fData.fTime.push_back( 0.01*fGtuLength*(Double_t)pix->GetGtu() ); // microseconds/100
fData.fHits.push_back( pix->GetCounts() );
fNumHits += pix->GetCounts();
}
fData.fNumPoints = fNumPoints;
fData.fNumHits = fNumHits;
fData.fGtuLength = fGtuLength;
Msg(EsafMsg::Info) << "Processing " << fNumPoints << " pixel and " << fNumHits << " hits" << MsgDispatch;
if(fData.fNumPoints < fNumPointsMin){
fQuality = -1;
Msg(EsafMsg::Warning) << "Not enough pixels selected ( minimum is "<< fNumPointsMin <<" ): terminate this event." << MsgDispatch;
return kFALSE;
} else {
fQuality = 0;
}
fTrueTheta = fEv->GetHeader().GetTrueTheta();
fTruePhi = fEv->GetHeader().GetTruePhi();
fTrueDir.SetMagThetaPhi( 1., fTrueTheta, fTruePhi);
Msg(EsafMsg::Info) << "TRUTH: Theta = " << fTrueTheta*DEG << " deg, Phi = " << fTruePhi*DEG << " deg, Energy = " << header.GetTrueEnergy() << " MeV." << MsgDispatch;
TVector3 fTrueMaxOldsdr = fEv->GetHeader().GetTrueShowerMaxPos();
Double_t xmax=fTrueMaxOldsdr.x()/km;
Double_t ymax=fTrueMaxOldsdr.y()/km;
Double_t zmax=fHISS-(fTrueMaxOldsdr.z()/km);
fTrueMax.SetXYZ(xmax,ymax,zmax);
Double_t phiMaxTmp=fTrueMax.Phi();
fTrueMax.SetPhi(phiMaxTmp+PI);
fTrueNorm = fTrueMax.Cross(fTrueDir);
if ( fTrueNorm.Z() > 0 ) fTrueNorm *= -1;
Bool_t fShapeSelectionDone = kFALSE;
if ( fDoHough ) {
UseHoughandFindPlane();
fShapeSelectionDone = kFALSE;
} else if ( !fDoShapeSelection ) {
FindPlane();
fShapeSelectionDone = kFALSE;
} else if ( fDoShapeSelection ) {
UseShapeandFindPlane();
if( fNumPointsSel < fNumPointsMin ){
Msg(EsafMsg::Info) << "Impossible to use selection by shape: only " << fNumPointsSel << " pixels selected!" <<MsgDispatch;
fShapeSelectionDone = kFALSE;
FindPlane();
} else { fShapeSelectionDone = kFALSE; }
}
Double_t fangleDirNorm = fTrueDir.Angle(fNorm);
fDeltaTDP = fTrueNorm.Angle(fNorm);
Msg(EsafMsg::Info) << "Found TrackDetectorPlane with error = " << fDeltaTDP*DEG << " deg, fangleDirNorm = " << fangleDirNorm*DEG << " deg" << MsgDispatch;
fW = (TVector3(0,0,1).Cross(fNorm)).Unit();
fU = (fNorm.Cross(fW)).Unit();
if( fOptionSelectionHough || (fUseShapeSelectioninModules && fShapeSelectionDone ) ){
fNumPoints = fNumPointsSel;
fNumHits = fNumHitsSel;
fData.fNumPoints = fNumPoints;
fData.fNumHits = fNumHits;
fData.fSpPos.clear();
fData.fTime.clear();
fData.fHits.clear();
fData.fSigmaTheta.clear();
fData.fSigmaPhi.clear();
for(Int_t i=0;i<fNumPointsSel;i++){
fData.fSpPos.push_back(fData.fSpPosSel[i]);
fData.fTime.push_back(fData.fTimeSel[i]);
fData.fHits.push_back(fData.fHitsSel[i]);
fData.fSigmaTheta.push_back(fData.fSigmaThetaSel[i]);
fData.fSigmaPhi.push_back(fData.fSigmaPhiSel[i]);
}
}
switch (fMethodIdentifier) {
case 1: AA1(); break;
case 2: NE1(); break;
case 3: AA2(); break;
case 4: NE2(); break;
case 5: AE1(); break;
case 6: all(); break;
default: Msg(EsafMsg::Panic) << "Wrong fMethod config value in TrackDirection2Module.cfg" << MsgDispatch; break;
}
fData.Clear();
//save data
MyData()->Add("NormX",fNorm.X());
MyData()->Add("NormY",fNorm.Y());
MyData()->Add("NormZ",fNorm.Z());
MyData()->Add("WAxisX",fW.X());
MyData()->Add("WAxisY",fW.Y());
MyData()->Add("WAxisZ",fW.Z());
MyData()->Add("Theta",fTHETAreco);
MyData()->Add("Phi",fPHIreco);
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::PostProcess() {
//
// Post-processing method
//
if( fQuality == 0) {
fRecoEventsCounter++;
vecDeltaTDP.push_back(fDeltaTDP*DEG);
if( fMethodIdentifier == 6 ) {
vecDeltaEASDirAA1.push_back(fDeltaEASDirAA1*DEG);
vecDeltaEASDirAA2.push_back(fDeltaEASDirAA2*DEG);
vecDeltaEASDirAE1.push_back(fDeltaEASDirAE1*DEG);
vecDeltaEASDirNE1.push_back(fDeltaEASDirNE1*DEG);
vecDeltaEASDirNE2.push_back(fDeltaEASDirNE2*DEG);
} else {
vecDeltaEASDir.push_back(fDeltaEASDir*DEG);
}
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::Done() {
//
// Module done method. Do some statistics about the renconstructed events
//
Msg(EsafMsg::Info) << "Number of reconstructed events = " << fRecoEventsCounter << MsgDispatch;
if (fRecoEventsCounter < 2) return kTRUE;
gDirectory->cd("/");
string pathname = "TD2Resolution";
TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
path->cd();
string aa1 = "AA1";
string aa2 = "AA2";
string ae1 = "AE1";
string ne1 = "NE1";
string ne2 = "NE2";
string tdp = "TDP";
DoStat(vecDeltaTDP,tdp);
if ( fMethodIdentifier == 6 ){
DoStat(vecDeltaEASDirAA1,aa1);
DoStat(vecDeltaEASDirAA2,aa2);
DoStat(vecDeltaEASDirAE1,ae1);
DoStat(vecDeltaEASDirNE1,ne1);
DoStat(vecDeltaEASDirNE2,ne2);
} else {
if ( fMethodIdentifier == 1 ) DoStat(vecDeltaEASDir,aa1);
if ( fMethodIdentifier == 2 ) DoStat(vecDeltaEASDir,ne1);
if ( fMethodIdentifier == 3 ) DoStat(vecDeltaEASDir,aa2);
if ( fMethodIdentifier == 4 ) DoStat(vecDeltaEASDir,ne2);
if ( fMethodIdentifier == 5 ) DoStat(vecDeltaEASDir,ae1);
}
vecDeltaEASDir.clear();
vecDeltaEASDirAA1.clear();
vecDeltaEASDirAA2.clear();
vecDeltaEASDirAE1.clear();
vecDeltaEASDirNE1.clear();
vecDeltaEASDirNE2.clear();
vecDeltaTDP.clear();
fRecoEventsCounter = 0;
Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
void TrackDirection2Module::DoStat(vector<Double_t> err, string namemethod){
//
// Statistics of reconstructed events.
//
Int_t best(0);
Int_t med(0);
Int_t worst(0);
Int_t numbins=400;
Int_t maxDeltaEAS=20;
TH1F *hR=new TH1F("hR","hR",numbins,0,maxDeltaEAS);
for( Int_t i=0; i<fRecoEventsCounter; i++ ) {
hR->Fill(err[i]);
if ( err[i] < fStat1 ) best++;
if ( err[i] >= fStat1 && err[i] <= fStat2 ) med++;
if ( err[i] > fStat2 ) worst++;
}
Bool_t foundresolution = kFALSE;
Stat_t under68(0);
Float_t frac(0);
Double_t resolution(0);
for(Int_t i=0;i<numbins;i++){
under68 += hR->GetBinContent(i);
frac=100*(Float_t)under68/((Float_t)fRecoEventsCounter);
if(frac>68){
resolution=i*maxDeltaEAS/((Float_t)numbins);
Msg(EsafMsg::Info)<<"Eas-Dir reconstruction ("<<namemethod.c_str()<<"):tResolution(68%)="<<resolution<<" deg"<<MsgDispatch;
foundresolution = kTRUE;
break;
}
}
if (!foundresolution)
Msg(EsafMsg::Info)<<"Eas-Dir reconstruction (" << namemethod.c_str() <<"):tResolution(68%) not found, too bad!!" << MsgDispatch;
string titleh = namemethod.c_str();
titleh += Form(" Resolution(68%)=%g deg", resolution);
hR->SetTitle(titleh.c_str());
hR->GetXaxis()->SetTitle("#Psi (deg)");
hR->Write();
delete hR;
Msg(EsafMsg::Info) << "tEvents with error under " << fStat1<< " deg = " <<best<< " ( " << ((Float_t)best/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
Msg(EsafMsg::Info) << "tEvents with error between "<<fStat1<<" and "<<fStat2<<" deg = " <<med<< " ( " << ((Float_t)med/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
Msg(EsafMsg::Info) << "tEvents with error above "<<fStat2<<" deg = " <<worst<< " ( " << ((Float_t)worst/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
return;
}
//_____________________________________________________________________________
void TrackDirection2Module::UserMemoryClean() {
//
// User memory clean
//
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::SaveRootData(RecoRootEvent *fRecoRootEvent) {
//
// Save data in the reco rootfile
//
fRecoRootEvent->GetRecoTrackDirection2().SetErrorTDP(fDeltaTDP*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetQuality(fQuality); // -1 if reconstruction is not possible in this context, otherwise 0
if(fMethodIdentifier == 6){
fRecoRootEvent->GetRecoTrackDirection2().SetAA1(fTHETArecoAA1*DEG, fPHIrecoAA1*DEG, fDeltaEASDirAA1*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetAA2(fTHETArecoAA2*DEG, fPHIrecoAA2*DEG, fDeltaEASDirAA2*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetNE1(fTHETArecoNE1*DEG, fPHIrecoNE1*DEG, fDeltaEASDirNE1*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetNE2(fTHETArecoNE2*DEG, fPHIrecoNE2*DEG, fDeltaEASDirNE2*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetAE1(fTHETArecoAE1*DEG, fPHIrecoAE1*DEG, fDeltaEASDirAE1*DEG);
}else{
fRecoRootEvent->GetRecoTrackDirection2().SetTheta(fTHETAreco*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetErrorTheta(fDeltaTheta*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetPhi(fPHIreco*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetErrorPhi(fDeltaPhi*DEG);
fRecoRootEvent->GetRecoTrackDirection2().SetErrorDirection(fDeltaEASDir*DEG);
}
return kTRUE;
}
//______________________________________________________________________________
void TrackDirection2Module::UseHoughandFindPlane(){
//
// Find track-detector plane (TDP) using Hough Transform
//
Msg(EsafMsg::Info) << "UseHoughandFindPlane()" << MsgDispatch;
vector<Double_t> x;
vector<Double_t> y;
vector<Double_t> t;
vector<Int_t> c;
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
t.push_back(fData.fTime[i]);
x.push_back(dummy.x()/dummy.z());
y.push_back(dummy.y()/dummy.z());
c.push_back(fData.fHits[i]);
}
// errors for the Hough fit
Float_t ex = 0.001;
Float_t ey = 0.001;
Float_t et = 2.5*0.01;
// use Hough fit to select track points
// x-t Hough fit
HoughFit *xtFit = new HoughFit( t, x, c, fOptionSelectionHough,et,ex,fPointsId);
Double_t vX = xtFit->GetSlope();
Double_t wX = xtFit->GetOffset();
vector<Int_t> idSelX = xtFit->GetIdSel();
delete xtFit;
// y-t Hough fit
HoughFit *ytFit = new HoughFit( t, y, c, fOptionSelectionHough,et,ey,fPointsId);
Double_t vY = ytFit->GetSlope();
Double_t wY = ytFit->GetOffset();
vector<Int_t> idSelY = ytFit->GetIdSel();
delete ytFit;
// compute the normal vector to the TDP
fNorm.SetXYZ(vY,-vX,wY*vX-wX*vY);
fNorm.SetMag(1.);
if ( fNorm.Z() > 0 ) fNorm *= -1;
//use the following only to debug (to see how the projection on plane Z=1 versus time appear)
if( fDoGraphUseShape ){
Double_t xprog[fNumPoints], yprog[fNumPoints], tprog[fNumPoints];
Double_t xline[fNumPoints],yline[fNumPoints];
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
tprog[i] = fData.fTime[i];
xline[i] = vX*fData.fTime[i] + wX;
yline[i]=vY*fData.fTime[i]+wY;
xprog[i] = dummy.x()/dummy.z();
yprog[i] = dummy.y()/dummy.z();
}
TGraph *gxypro = new TGraph(fNumPoints,yprog,xprog);
gxypro->SetNameTitle("gxypro","xprog vs yprog");gxypro->SetMarkerSize(.5);gxypro->SetMarkerStyle(23);
TGraph *gxpro = new TGraph(fNumPoints,tprog,xprog);
gxpro->SetNameTitle("gxpro","xprog vs time");gxpro->SetMarkerSize(.5);gxpro->SetMarkerStyle(23);
TGraph *gypro = new TGraph(fNumPoints,tprog,yprog);
gypro->SetNameTitle("gypro","yprog vs time");gypro->SetMarkerSize(.5);gypro->SetMarkerStyle(23);
TGraph *gxproretta = new TGraph(fNumPoints,tprog,xline);
gxproretta->SetMarkerColor(4);gxproretta->SetLineColor(4);
TGraph *gyproretta = new TGraph(fNumPoints,tprog,yline);
gyproretta->SetMarkerColor(4);gyproretta->SetLineColor(4);
TCanvas *cUSFP=new TCanvas("cUSFP","cUSFP",1);
cUSFP->Divide(3,1);
cUSFP->cd(1);gxpro->Draw("AP");gxproretta->Draw("PLsame");
cUSFP->cd(2);gypro->Draw("AP");gyproretta->Draw("PLsame");
cUSFP->cd(3);gxypro->Draw("AP");
cUSFP->Write();
}
//cout<<"fPointsId.size()="<<fPointsId.size()<<"tidSelX.size()="<<idSelX.size()<<"tidSelY.size()="<<idSelY.size()<<endl;
vector<Int_t> idSelectedMerged(idSelY.size()+idSelX.size());
merge(idSelX.begin(), idSelX.end(),idSelY.begin(), idSelY.end(),idSelectedMerged.begin());
vector<Int_t>::iterator different = unique(idSelectedMerged.begin(),idSelectedMerged.end());
idSelectedMerged.erase(different, idSelectedMerged.end());
fNumPointsSel = idSelectedMerged.size();
Msg(EsafMsg::Info) << "Selected " << fNumPointsSel << " points over " << fPointsId.size() << MsgDispatch;
// if there are > 10 points selected recompute the TDP more precisely
fNumHitsSel = 0;
if( fNumPointsSel > 10 ){
vector<Int_t> csel;
vector<Double_t> xsel;
vector<Double_t> ysel;
vector<Double_t> errxsel;
vector<Double_t> errysel;
vector<Double_t> tsel;
vector<Double_t> xsel1;
vector<Double_t> ysel1;
vector<Double_t> tsel1;
for(Int_t i=0; i<fNumPointsSel; i++) {
RecoPixelData *pix = fEv->GetRecoPixelData( idSelectedMerged[i] );
TVector3 pos(0,0,0);
Double_t phitmp = pix->GetPhi() + PI;
pos.SetMagThetaPhi( 1, pix->GetTheta(), phitmp );
xsel.push_back( pos.x()/pos.z() );
ysel.push_back( pos.y()/pos.z() );
tsel.push_back( 0.01*(Double_t)2.500*pix->GetGtu() );
csel.push_back( pix->GetCounts() );
fNumHitsSel += csel[i];
Double_t errphi = pix->GetPhiSigma();
Double_t errtheta = pix->GetThetaSigma();
fData.fSpPosSel.push_back(pos);
fData.fTimeSel.push_back(tsel[i]);
fData.fHitsSel.push_back(csel[i]);
fData.fSigmaThetaSel.push_back(errtheta);
fData.fSigmaPhiSel.push_back(errphi);
Double_t dx = DeltaX(pos.Theta(), pos.Phi(), errtheta, errphi);
Double_t dy = DeltaY(pos.Theta(), pos.Phi(), errtheta, errphi);
errxsel.push_back(dx/csel[i]);
errysel.push_back(dy/csel[i]);
// add hits for median fit
for (Int_t j(0); j<pix->GetCounts(); j++) {
xsel1.push_back( pos.x()/pos.z() );
ysel1.push_back( pos.y()/pos.z() );
tsel1.push_back( 0.01*(Double_t)2.500*pix->GetGtu() );
}
}
TVector3 fNormA, fNormB, fNormC;
Double_t vX2(0), vY2(0), wX2(0), wY2(0);
if ( fDebugInfo || useHFplane ) { // recalculate normal to TDP using Hough fit
HoughFit *xhf1 = new HoughFit( tsel, xsel, csel, 0,et,ex);
vX2 = xhf1->GetSlope();
wX2 = xhf1->GetOffset();
delete xhf1;
HoughFit *yhf1 = new HoughFit( tsel, ysel, csel, 0,et,ey);
vY2 = yhf1->GetSlope();
wY2 = yhf1->GetOffset();
delete yhf1;
fNormA.SetXYZ(vY2,-vX2,wY2*vX2-wX2*vY2);
fNormA.SetMag(1.);
if ( fNormA.Z() > 0 ) fNormA *= -1;
Double_t fDeltaTDPhh = fTrueNorm.Angle(fNormA);
Msg(EsafMsg::Info) << "TDP error (hough selection + hough fit) = " << fDeltaTDPhh*DEG << MsgDispatch;
}
if ( fDebugInfo || useLFplane ) { // recalculate normal to TDP using least squares fit
LeastSquaresFit *xlf1 = new LeastSquaresFit( fNumPointsSel, tsel, xsel, errxsel );
Double_t vX2ls = xlf1->GetSlope();
Double_t wX2ls = xlf1->GetOffset();
delete xlf1;
LeastSquaresFit *ylf1 = new LeastSquaresFit( fNumPointsSel, tsel, ysel, errysel );
Double_t vY2ls = ylf1->GetSlope();
Double_t wY2ls = ylf1->GetOffset();
delete ylf1;
fNormB.SetXYZ(vY2ls,-vX2ls,wY2ls*vX2ls-wX2ls*vY2ls);
fNormB.SetMag(1.); if ( fNormB.Z() > 0 ) fNormB *= -1;
Double_t fDeltaTDPhls = fTrueNorm.Angle(fNormB);
Msg(EsafMsg::Info) << "TDP error (hough selection + linear fit) = " << fDeltaTDPhls*DEG << MsgDispatch;
}
if ( fDebugInfo || useHFplane ) { // recalculate normal to TDP using median fit
MedianFit *xmf1 = new MedianFit( fNumHitsSel, tsel1, xsel1 ); // CHANGED
Double_t vX2m = xmf1->GetSlope();
Double_t wX2m = xmf1->GetOffset();
delete xmf1;
MedianFit *ymf1 = new MedianFit( fNumHitsSel, tsel1, ysel1 ); // CHANGED
Double_t vY2m = ymf1->GetSlope();
Double_t wY2m = ymf1->GetOffset();
delete ymf1;
fNormC.SetXYZ(vY2m,-vX2m,wY2m*vX2m-wX2m*vY2m);
fNormC.SetMag(1.);
if ( fNormC.Z() > 0 ) fNormC *= -1;
Double_t fDeltaTDPhm = fTrueNorm.Angle(fNormC);
Msg(EsafMsg::Info) << "TDP error (hough selection + median fit) = " << fDeltaTDPhm*DEG << MsgDispatch;
}
if ( useLFplane ) fNorm = fNormB;
else if ( useMFplane ) fNorm = fNormC;
else if ( useHFplane ) fNorm = fNormA;
fNorm.SetMag(1.);
if ( fNorm.Z() > 0 ) fNorm *= -1;
//use the following only to debug (to see how the projection on plane Z=1 versus time appear)
if( fDoGraphUseShape ){
Double_t xprog2[fNumPointsSel],yprog2[fNumPointsSel],tprog2[fNumPointsSel],xline2[fNumPointsSel],yline2[fNumPointsSel];
for( Int_t i=0; i<fNumPointsSel; i++ ) {
tprog2[i] = tsel[i];
xline2[i] = vX2*tsel[i]+wX2;
yline2[i] = vY2*tsel[i]+wY2;
xprog2[i] = xsel[i];
yprog2[i] = ysel[i];
}
TGraph *gxypro2 = new TGraph(fNumPointsSel,yprog2,xprog2);
gxypro2->SetNameTitle("gxypro2","xprog2 vs yprog2");gxypro2->SetMarkerSize(.5);gxypro2->SetMarkerStyle(23);
TGraph *gxpro2 = new TGraph(fNumPointsSel,tprog2,xprog2);
gxpro2->SetNameTitle("gxpro2","xprog2 vs time2");gxpro2->SetMarkerSize(.5);gxpro2->SetMarkerStyle(23);
TGraph *gypro2 = new TGraph(fNumPointsSel,tprog2,yprog2);
gypro2->SetNameTitle("gypro2","yprog2 vs time2");gypro2->SetMarkerSize(.5);gypro2->SetMarkerStyle(23);
TGraph *gxproretta2 = new TGraph(fNumPointsSel,tprog2,xline2);
gxproretta2->SetMarkerColor(4);gxproretta2->SetLineColor(4);
TGraph *gyproretta2 = new TGraph(fNumPointsSel,tprog2,yline2);
gyproretta2->SetMarkerColor(4);gyproretta2->SetLineColor(4);
TCanvas *cUSFP2=new TCanvas("cUSFP2","cUSFP2",1);
cUSFP2->Divide(3,1);
cUSFP2->cd(1);gxpro2->Draw("AP");gxproretta2->Draw("PLsame");
cUSFP2->cd(2);gypro2->Draw("AP");gyproretta2->Draw("PLsame");
cUSFP2->cd(3);gxypro2->Draw("AP");
cUSFP2->Write();
}
}
}
//______________________________________________________________________________
void TrackDirection2Module::FindPlane(){
//
// Find the track-detector plane (TDP)
//
Msg(EsafMsg::Info) << "FindPlane()" << MsgDispatch;
vector<Double_t> xpro;
vector<Double_t> ypro;
vector<Double_t> errxpro;
vector<Double_t> errypro;
vector<Double_t> tpro;
vector<Int_t> cpro;
vector<Double_t> xpro1;
vector<Double_t> ypro1;
vector<Double_t> tpro1;
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
xpro.push_back(dummy.X()/dummy.Z());
ypro.push_back(dummy.Y()/dummy.Z());
cpro.push_back(fData.fHits[i]);
tpro.push_back(fData.fTime[i]);
errxpro.push_back(dx/fData.fHits[i]);
errypro.push_back(dy/fData.fHits[i]);
// add hits for median fit
for( Int_t j(0); j<fData.fHits[i]; j++ ) {
xpro1.push_back(dummy.X()/dummy.Z());
ypro1.push_back(dummy.Y()/dummy.Z());
tpro1.push_back(fData.fTime[i]);
}
}
TVector3 fNormA, fNormB, fNormC;
if ( fDebugInfo || useHFplane ) { // calculate normal to TDP using Hough fit
// errors for the Hough fit
Float_t ex = 0.001;
Float_t ey = 0.001;
Float_t et = 2.5*0.01;
HoughFit *xhf1 = new HoughFit( tpro, xpro, cpro, 0,et,ex);
Double_t vX2 = xhf1->GetSlope();
Double_t wX2 = xhf1->GetOffset();
delete xhf1;
HoughFit *yhf1 = new HoughFit( tpro, ypro, cpro, 0,et,ey);
Double_t vY2 = yhf1->GetSlope();
Double_t wY2 = yhf1->GetOffset();
delete yhf1;
fNormA.SetXYZ(vY2,-vX2,wY2*vX2-wX2*vY2);
fNormA.SetMag(1.);
if ( fNormA.Z() > 0 ) fNormA *= -1;
Double_t fDeltaTDPhh = fTrueNorm.Angle(fNormA);
Msg(EsafMsg::Info) << "TDP error (hough fit) = " << fDeltaTDPhh*DEG << MsgDispatch;
}
if ( fDebugInfo || useLFplane ) { // calculate normal to TDP using least squares fit
LeastSquaresFit *xlf1 = new LeastSquaresFit( fNumPoints, tpro, xpro, errxpro );
Double_t vX2ls = xlf1->GetSlope();
Double_t wX2ls = xlf1->GetOffset();
delete xlf1;
LeastSquaresFit *ylf1 = new LeastSquaresFit( fNumPoints, tpro, ypro, errypro );
Double_t vY2ls = ylf1->GetSlope();
Double_t wY2ls = ylf1->GetOffset();
delete ylf1;
fNormB.SetXYZ(vY2ls,-vX2ls,wY2ls*vX2ls-wX2ls*vY2ls);
fNormB.SetMag(1.); if ( fNormB.Z() > 0 ) fNormB *= -1;
Double_t fDeltaTDPhls = fTrueNorm.Angle(fNormB);
Msg(EsafMsg::Info) << "TDP error (linear fit) = " << fDeltaTDPhls*DEG << MsgDispatch;
}
if ( fDebugInfo || useHFplane ) { // recalculate normal to TDP using median fit
MedianFit *xmf1 = new MedianFit( fNumHits, tpro1, xpro1 );
Double_t vX2m = xmf1->GetSlope();
Double_t wX2m = xmf1->GetOffset();
delete xmf1;
MedianFit *ymf1 = new MedianFit( fNumHits, tpro1, ypro1 );
Double_t vY2m = ymf1->GetSlope();
Double_t wY2m = ymf1->GetOffset();
delete ymf1;
fNormC.SetXYZ(vY2m,-vX2m,wY2m*vX2m-wX2m*vY2m);
fNormC.SetMag(1.);
if ( fNormC.Z() > 0 ) fNormC *= -1;
Double_t fDeltaTDPhm = fTrueNorm.Angle(fNormC);
Msg(EsafMsg::Info) << "TDP error (median fit) = " << fDeltaTDPhm*DEG << MsgDispatch;
}
if ( useLFplane ) fNorm = fNormB;
else if ( useMFplane ) fNorm = fNormC;
else if ( useHFplane ) fNorm = fNormA;
fNorm.SetMag(1.);
if ( fNorm.Z() > 0 ) fNorm *= -1;
cpro.clear();
xpro.clear();
ypro.clear();
tpro.clear();
errxpro.clear();
errypro.clear();
}
//______________________________________________________________________________
void TrackDirection2Module::UseShapeandFindPlane(){
//
// Find the track-detector plane (TDP) using the shape method
//
vector<Int_t> cpro;
vector<Double_t> xpro;
vector<Double_t> ypro;
vector<Double_t> errxpro;
vector<Double_t> errypro;
vector<Double_t> tpro;
vector<Double_t> xpro1;
vector<Double_t> ypro1;
vector<Double_t> tpro1;
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
xpro.push_back(dummy.X()/dummy.Z());
ypro.push_back(dummy.Y()/dummy.Z());
tpro.push_back(fData.fTime[i]);
cpro.push_back(fData.fHits[i]);
errxpro.push_back(dx/fData.fHits[i]);
errypro.push_back(dy/fData.fHits[i]);
// add hits for median fit
for(Int_t j(0); j<fData.fHits[i]; i++) {
xpro1.push_back(dummy.X()/dummy.Z());
ypro1.push_back(dummy.Y()/dummy.Z());
tpro1.push_back(fData.fTime[i]);
}
}
// FIRST STEP: median fit of points
MedianFit *xMf1 = new MedianFit(fNumHits, tpro1, xpro1); //CHANGED
Double_t sX = xMf1->GetSlope();
Double_t oX = xMf1->GetOffset();
Double_t aX = xMf1->GetAbsoluteDeviation();
delete xMf1;
MedianFit *yMf1 = new MedianFit(fNumHits, tpro1, ypro1); //CHANGED
Double_t sY = yMf1->GetSlope();
Double_t oY = yMf1->GetOffset();
Double_t aY = yMf1->GetAbsoluteDeviation();
delete yMf1;
//use the following only to debug (to see how the projection on plane Z=1 versus time appear)
if( fDoGraphUseShape ){
Double_t xprog[fNumPoints], yprog[fNumPoints], tprog[fNumPoints];
Double_t xline[fNumPoints], yline[fNumPoints];
Double_t xlineup[fNumPoints], ylineup[fNumPoints];
Double_t xlinedown[fNumPoints], ylinedown[fNumPoints];
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
tprog[i] = fData.fTime[i];
xline[i] = sX * fData.fTime[i] + oX;
xlineup[i] = xline[i] + fMulti1 * aX;
xlinedown[i] = xline[i] - fMulti1 * aX;
yline[i] = sY * fData.fTime[i] + oY;
ylineup[i] = yline[i] + fMulti1 * aY;
ylinedown[i] = yline[i] - fMulti1 * aY;
xprog[i] = dummy.X()/dummy.Z();
yprog[i] = dummy.Y()/dummy.Z();
}
TGraph *gxypro = new TGraph(fNumPoints,yprog,xprog);
gxypro->SetNameTitle("gxypro","xprog vs yprog");gxypro->SetMarkerSize(.5);gxypro->SetMarkerStyle(23);
TGraph *gxpro = new TGraph(fNumPoints,tprog,xprog);
gxpro->SetNameTitle("gxpro","xprog vs time");gxpro->SetMarkerSize(.5);gxpro->SetMarkerStyle(23);
TGraph *gypro = new TGraph(fNumPoints,tprog,yprog);
gypro->SetNameTitle("gypro","yprog vs time");
gypro->SetMarkerSize(.5);
gypro->SetMarkerStyle(23);
TGraph *gxproretta = new TGraph(fNumPoints,tprog,xline);
gxproretta->SetMarkerColor(4);gxproretta->SetLineColor(4);
TGraph *gyproretta = new TGraph(fNumPoints,tprog,yline);
gyproretta->SetMarkerColor(4);gyproretta->SetLineColor(4);
TGraph *gxprorettaup = new TGraph(fNumPoints,tprog,xlineup);
gxprorettaup->SetMarkerColor(2);gxprorettaup->SetLineColor(2);
TGraph *gyprorettaup = new TGraph(fNumPoints,tprog,ylineup);
gyprorettaup->SetMarkerColor(2);gyprorettaup->SetLineColor(2);
TGraph *gxprorettadown = new TGraph(fNumPoints,tprog,xlinedown);
gxprorettadown->SetMarkerColor(2);gxprorettadown->SetLineColor(2);
TGraph *gyprorettadown = new TGraph(fNumPoints,tprog,ylinedown);
gyprorettadown->SetMarkerColor(2);gyprorettadown->SetLineColor(2);
TCanvas *cUSFP=new TCanvas("cUSFP","cUSFP",1);
cUSFP->Divide(3,1);
cUSFP->cd(1);gxpro->Draw("AP");gxproretta->Draw("PLsame");gxprorettaup->Draw("PLsame");gxprorettadown->Draw("PLsame");
cUSFP->cd(2);gypro->Draw("AP");gyproretta->Draw("PLsame");gyprorettaup->Draw("PLsame");gyprorettadown->Draw("PLsame");
cUSFP->cd(3);gxypro->Draw("AP");
cUSFP->Write();
}
// make a shape selection of track points
vector<Double_t> xprosel;
vector<Double_t> yprosel;
vector<Double_t> errxprosel;
vector<Double_t> erryprosel;
vector<Double_t> tprosel;
vector<Double_t> xprosel1;
vector<Double_t> yprosel1;
vector<Double_t> tprosel1;
vector<TVector3> fSpPosSeltmp;
vector<Double_t> fTimeSeltmp;
vector<Int_t> fHitsSeltmp;
vector<Double_t> fSigmaThetaSeltmp;
vector<Double_t> fSigmaPhiSeltmp;
Int_t fNumPointsSeltmp = 0;
Int_t fNumHitsSeltmp = 0;
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
Double_t xprotmp=dummy.X()/dummy.Z();
Double_t yprotmp=dummy.Y()/dummy.Z();
Double_t tprotmp=fData.fTime[i];
// shape selection with multiplier fMulti1
if( Abs(xprotmp-sX*tprotmp-oX) < fMulti1*aX || Abs(yprotmp-sY*tprotmp-oY) < fMulti1*aY) {
Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
fSpPosSeltmp.push_back(dummy);
fTimeSeltmp.push_back(fData.fTime[i]);
fHitsSeltmp.push_back(fData.fHits[i]);
fSigmaThetaSeltmp.push_back(fData.fSigmaTheta[i]);
fSigmaPhiSeltmp.push_back(fData.fSigmaPhi[i]);
xprosel.push_back(xprotmp);
yprosel.push_back(yprotmp);
tprosel.push_back(tprotmp);
errxprosel.push_back(dx/fData.fHits[i]);
erryprosel.push_back(dy/fData.fHits[i]);
fNumPointsSeltmp++;
fNumHitsSeltmp += fData.fHits[i];
//add hits for median fit
for( Int_t j(0); j<fData.fHits[i]; j++) {
xprosel1.push_back(xprotmp);
yprosel1.push_back(yprotmp);
tprosel1.push_back(tprotmp);
}
}
}
Msg(EsafMsg::Info) << "UseShapeandFindPlane(), first step: selected " << fNumPointsSeltmp << " from " << fNumPoints << " pixels" << MsgDispatch;
if ( fNumPointsSeltmp < fNumPointsMin ){ // return if not enough points
fNumPointsSel = fNumPointsSeltmp;
return;
}
// SECOND STEP: new median fit of selected points
MedianFit *xMf2 = new MedianFit(fNumHitsSeltmp, tprosel1, xprosel1); //CHANGED
Double_t sX2 = xMf2->GetSlope();
Double_t oX2 = xMf2->GetOffset();
Double_t aX2 = xMf2->GetAbsoluteDeviation();
delete xMf2;
MedianFit *yMf2 = new MedianFit(fNumHitsSeltmp, tprosel1, yprosel1); //CHANGED
Double_t sY2 = yMf2->GetSlope();
Double_t oY2 = yMf2->GetOffset();
Double_t aY2 = yMf2->GetAbsoluteDeviation();
delete yMf2;
// new shape selection of points from median fit results
vector<Double_t> xprosel2;
vector<Double_t> yprosel2;
vector<Double_t> errxprosel2;
vector<Double_t> erryprosel2;
vector<Double_t> tprosel2;
vector<Int_t> cprosel2;
vector<Double_t> xprosel3;
vector<Double_t> yprosel3;
vector<Double_t> tprosel3;
fNumPointsSel = 0;
fNumHitsSel = 0;
for( Int_t i=0; i<fNumPointsSeltmp; i++ ) {
TVector3 dummy = fSpPosSeltmp[i];
Double_t xprotmp = dummy.X()/dummy.Z();
Double_t yprotmp = dummy.Y()/dummy.Z();
Double_t tprotmp = fTimeSeltmp[i];
// shape selction with multiplier fMulti2
if( Abs(xprotmp-sX2*tprotmp-oX2) < fMulti2*aX2 || TMath::Abs(yprotmp-sY2*tprotmp-oY2) < fMulti2*aY2 ){
Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
fData.fSpPosSel.push_back(dummy);
fData.fTimeSel.push_back(fTimeSeltmp[i]);
fData.fHitsSel.push_back(fHitsSeltmp[i]);
fData.fSigmaThetaSel.push_back(fSigmaThetaSeltmp[i]);
fData.fSigmaPhiSel.push_back(fSigmaPhiSeltmp[i]);
cprosel2.push_back(fHitsSeltmp[i]);
xprosel2.push_back(xprotmp);
yprosel2.push_back(yprotmp);
tprosel2.push_back(tprotmp);
errxprosel2.push_back(dx/fHitsSeltmp[i]);
erryprosel2.push_back(dy/fHitsSeltmp[i]);
fNumPointsSel++;
fNumHitsSel += fHitsSeltmp[i];
//add hits for median fit
for (Int_t j(0); j<fHitsSeltmp[i]; j++) {
xprosel3.push_back(xprotmp);
yprosel3.push_back(yprotmp);
tprosel3.push_back(tprotmp);
}
}
}
Msg(EsafMsg::Info) << "UseShapeandFindPlane(), second step: selected " << fNumPointsSel << " from " << fNumPointsSeltmp << " pixels" << MsgDispatch;
if ( fNumPointsSel < fNumPointsMin ) return; // return if not enough points
Double_t vX(0), wX(0), vY(0), wY(0);
// fit of selected points
if ( useLFplane ) {
LeastSquaresFit *xtFit = new LeastSquaresFit(fNumPointsSel, tprosel2, xprosel2, errxprosel2);
vX = xtFit->GetSlope();
wX = xtFit->GetOffset();
delete xtFit;
LeastSquaresFit *ytFit = new LeastSquaresFit(fNumPointsSel, tprosel2, yprosel2, erryprosel2);
vY = ytFit->GetSlope();
wY = ytFit->GetOffset();
delete ytFit;
}
if ( useMFplane ) {
MedianFit *xtFit = new MedianFit(fNumHitsSel, tprosel3, xprosel3);
vX = xtFit->GetSlope();
wX = xtFit->GetOffset();
delete xtFit;
MedianFit *ytFit = new MedianFit(fNumHitsSel, tprosel3, yprosel3);
vY = ytFit->GetSlope();
wY = ytFit->GetOffset();
delete ytFit;
}
if ( useHFplane ) {
Float_t ex = 0.001;
Float_t ey = 0.001;
Float_t et = 2.5*0.01;
HoughFit *xhf1 = new HoughFit( tpro, xpro, cpro, 0,et,ex);
vX = xhf1->GetSlope();
wX = xhf1->GetOffset();
delete xhf1;
HoughFit *yhf1 = new HoughFit( tpro, ypro, cpro, 0,et,ey);
vY = yhf1->GetSlope();
wY = yhf1->GetOffset();
delete yhf1;
}
// calculate normal versor to the TDP
fNorm.SetXYZ(vY,-vX,wY*vX-wX*vY);
fNorm.SetMag(1.);
if ( fNorm.Z() > 0 ) fNorm *= -1;
//use the following only to debug (to see how the projection on plane Z=1 versus time appear after selection)
if( fDoGraphUseShape ) {
Double_t xprog2[fNumPointsSel],yprog2[fNumPointsSel],tprog2[fNumPointsSel];
Double_t xline2[fNumPointsSel],yline2[fNumPointsSel];
//Double_t xline2up[fNumPoints],yline2up[fNumPoints];
//Double_t xline2down[fNumPoints],yline2down[fNumPoints];
for( Int_t i=0; i<fNumPointsSel; i++ ) {
TVector3 dummy = fData.fSpPosSel[i];
//xline2[i] = vX * fData.fTimeSel[i] + wX;
//yline2[i] = vY * fData.fTimeSel[i] + wY;
xline2[i] = vX * tprosel2[i] + wX;
yline2[i] = vY * tprosel2[i] + wY;
//xline2up[i] = xline2[i] + fMulti2 * a2X;
//xline2down[i] = xline2[i] - fMulti2 * a2X;
//yline2up[i] = yline2[i] + fMulti2 * a2Y;
//yline2down[i] = yline2[i] - fMulti2 * a2Y;
xprog2[i] = xprosel2[i]; // xprog2[i] = dummy.X()/dummy.Z();
yprog2[i] = yprosel2[i]; // yprog2[i] = dummy.Y()/dummy.Z();
tprog2[i] = tprosel2[i]; // tprog2[i] = fData.fTimeSel[i];
}
TGraph *gxypro2 = new TGraph(fNumPointsSel,yprog2,xprog2);
gxypro2->SetNameTitle("gxypro2","xprog2 vs yprog2");
gxypro2->SetMarkerSize(.5);
gxypro2->SetMarkerStyle(23);
TGraph *gxpro2 = new TGraph(fNumPointsSel,tprog2,xprog2);
gxpro2->SetNameTitle("gxpro2","xprog2 vs time2");
gxpro2->SetMarkerSize(.5);
gxpro2->SetMarkerStyle(23);
TGraph *gypro2 = new TGraph(fNumPointsSel,tprog2,yprog2);
gypro2->SetNameTitle("gypro2","yprog2 vs time2");
gypro2->SetMarkerSize(.5);
gypro2->SetMarkerStyle(23);
TGraph *gxproretta2 = new TGraph(fNumPointsSel,tprog2,xline2);
gxproretta2->SetMarkerColor(4);gxproretta2->SetLineColor(4);
TGraph *gyproretta2 = new TGraph(fNumPointsSel,tprog2,yline2);
gyproretta2->SetMarkerColor(4);gyproretta2->SetLineColor(4);
/*TGraph *gxproretta2up = new TGraph(fNumPointsSel,tprog2,xline2up);
gxproretta2up->SetMarkerColor(2);gxproretta2up->SetLineColor(2);
TGraph *gyproretta2up = new TGraph(fNumPointsSel,tprog2,yline2up);
gyproretta2up->SetMarkerColor(2);gyproretta2up->SetLineColor(2);
TGraph *gxproretta2down = new TGraph(fNumPointsSel,tprog2,xline2down);
gxproretta2down->SetMarkerColor(2);gxproretta2down->SetLineColor(2);
TGraph *gyproretta2down = new TGraph(fNumPointsSel,tprog2,yline2down);
gyproretta2down->SetMarkerColor(2);gyproretta2down->SetLineColor(2);*/
TCanvas *cUSFP2=new TCanvas("cUSFP2","cUSFP2",1);
cUSFP2->Divide(3,1);
cUSFP2->cd(1);gxpro2->Draw("AP");gxproretta2->Draw("PLsame");//gxproretta2up->Draw("PLsame");gxproretta2down->Draw("PLsame");
cUSFP2->cd(2);gypro2->Draw("AP");gyproretta2->Draw("PLsame");//gyproretta2up->Draw("PLsame");gyproretta2down->Draw("PLsame");
cUSFP2->cd(3);gxypro2->Draw("AP");
cUSFP2->Write();
}
// clear of containers
cpro.clear(); xpro.clear(); ypro.clear(); tpro.clear();
errxpro.clear(); errypro.clear();
xprosel.clear(); yprosel.clear(); tprosel.clear();
errxprosel.clear(); erryprosel.clear();
cprosel2.clear(); xprosel2.clear(); yprosel2.clear();
errxprosel2.clear(); erryprosel2.clear(); tprosel2.clear();
Msg(EsafMsg::Info) << "UseShapeandFindPlane(): " << fNumPointsSel << " pixel and " << fNumHitsSel << " hits selected " << MsgDispatch;
return;
}
//______________________________________________________________________________
void TrackDirection2Module::AA1() {
//
// ANALITIC APPROXIMATED 1 method.
// Shower constant angular velocity of the shower approximation.
// This method initialize the parameters for other all fit methods.
//
Msg(EsafMsg::Info) << "Using method AA1()" << MsgDispatch;
fData.fMedDir.SetXYZ(0,0,0);
fData.fMedTime = 0;
fData.fMedAlpha = 0;
vector<Double_t> erralpha;
vector<Int_t> counts;
//vectors of data for median fit
vector<Double_t> alphahits;
vector<Double_t> timehits;
Int_t numhits = 0;
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
Double_t alphatmp = ACos(dummy*fW);
fData.fAlpha.push_back(alphatmp);
erralpha.push_back(fErrAngle/Sqrt((Double_t)fData.fHits[i]));
counts.push_back(fData.fHits[i]);
numhits += fData.fHits[i];
fData.fMedTime += fData.fTime[i]*fData.fHits[i];
fData.fMedDir += fData.fSpPos[i]*fData.fHits[i];
//add hits for median fit
for( Int_t j(0); j<fData.fHits[i]; j++ ) {
alphahits.push_back(alphatmp);
timehits.push_back(fData.fTime[i]);
}
}
fData.fApparentTimeLenght = (Int_t)((fData.fTime[fNumPoints-1]-fData.fTime[0])/fGtuLength);
fData.fMedTime = fData.fMedTime/fNumHits;
fData.fMedDir = (fData.fMedDir*(1/(Double_t) fNumHits)).Unit();
fData.fMedAlpha = ACos(fData.fMedDir*fW);
// shower angular velocity fit: least squares, median and hough
Double_t fAngularSpeedLF(0), fQAngularSpeedLF(0), fAngularSpeedtodrawLF(0);
Double_t fAngularSpeedMF(0), fQAngularSpeedMF(0), fAngularSpeedtodrawMF(0);
Double_t fAngularSpeedHF(0), fQAngularSpeedHF(0), fAngularSpeedtodrawHF(0);
if ( fDebugInfo || fDoGraphUseShape || useLFaa1 ) {
LeastSquaresFit *ASl = new LeastSquaresFit(fNumPoints, fData.fTime, fData.fAlpha, erralpha);
fAngularSpeedLF = ASl->GetSlope()*0.01; // uom fix with C()
fQAngularSpeedLF = ASl->GetOffset();
fAngularSpeedtodrawLF = ASl->GetSlope();
delete ASl;
}
if ( fDebugInfo || fDoGraphUseShape || useMFaa1 ) {
MedianFit *ASm = new MedianFit(numhits, timehits, alphahits); //CHANGED
fAngularSpeedMF = ASm->GetSlope()*0.01; // uom fix with C()
fQAngularSpeedMF = ASm->GetOffset();
fAngularSpeedtodrawMF = ASm->GetSlope();
delete ASm;
}
if ( fDebugInfo || fDoGraphUseShape || useHFaa1 ) {
Float_t et = 0.01; // errors for the hough fit
Float_t ea = fErrAngle;
HoughFit *ASh = new HoughFit( fData.fTime, fData.fAlpha, counts, 0,et,ea);
fAngularSpeedHF = ASh->GetSlope()*0.01; // uom fix with C()
fQAngularSpeedHF = ASh->GetOffset();
fAngularSpeedtodrawHF = ASh->GetSlope();
delete ASh;
}
if (fDebugInfo)
Msg(EsafMsg::Info) << "method AA1(): fAngularSpeedLF = " << fAngularSpeedLF << " fAngularSpeedMF = " << fAngularSpeedMF << " fAngularSpeedHF = " << fAngularSpeedHF << MsgDispatch;
// use the following only to debug (plot angles alphai vs time)
Double_t alphai[fNumPoints],timei[fNumPoints],alphailinelf[fNumPoints],alphailinemf[fNumPoints],alphailinehf[fNumPoints];
if( fDoGraphUseShape ){
for( Int_t i=0; i<fNumPoints; i++ ) {
alphai[i]=fData.fAlpha[i]*DEG;
timei[i]=fData.fTime[i];
alphailinelf[i]=(fData.fTime[i]*fAngularSpeedtodrawLF+fQAngularSpeedLF)*DEG;
alphailinemf[i]=(fData.fTime[i]*fAngularSpeedtodrawMF+fQAngularSpeedMF)*DEG;
alphailinehf[i]=(fData.fTime[i]*fAngularSpeedtodrawHF+fQAngularSpeedHF)*DEG;
}
TGraph *gang = new TGraph(fNumPoints,timei,alphai);
gang->SetNameTitle("gang","alphai vs time;time;alpha");
gang->SetMarkerSize(.5);gang->SetMarkerStyle(23);
TGraph *ganglinelf = new TGraph(fNumPoints,timei,alphailinelf);ganglinelf->SetLineColor(4);
TGraph *ganglinemf = new TGraph(fNumPoints,timei,alphailinemf);ganglinemf->SetLineColor(6);
TGraph *ganglinehf = new TGraph(fNumPoints,timei,alphailinehf);ganglinehf->SetLineColor(8);
TCanvas *cAA1=new TCanvas("cAA1","cAA1",1);cAA1->Divide(1,1);
cAA1->cd(1);
gang->Draw("AP");
ganglinelf->Draw("PLsame");
ganglinemf->Draw("PLsame");
ganglinehf->Draw("PLsame");
TLegend *legend = new TLegend(0.8,0.9,1,0.6);
legend->AddEntry(ganglinelf,"least squares fit","l");
legend->AddEntry(ganglinemf,"median fit","l");
legend->AddEntry(ganglinehf,"hough fit","l");
legend->Draw();
cAA1->Write();
delete cAA1;
delete gang;
delete ganglinelf;
delete ganglinemf;
delete ganglinehf;
delete legend;
}
// selection of fit method
if ( useLFaa1 ) {
fAngularSpeed = fAngularSpeedLF;
} else if ( useMFaa1 ) {
fAngularSpeed = fAngularSpeedMF;
} else if ( useHFaa1 ) {
fAngularSpeed = fAngularSpeedHF;
} else {
fAngularSpeed = fAngularSpeedLF;
Msg(EsafMsg::Warning) << "AA1() fith method non specified in config file. Saving least squares fit results." << MsgDispatch;
}
// check if the result of the fit is a 'not a number'
if ( IsNaN(fAngularSpeed) ) {
Msg(EsafMsg::Warning) << "NaN angular speed, skipping this event" << MsgDispatch;
fAA1done = kTRUE;
fAA1nan = kTRUE;
return;
}
fHmax = 5; // HMax initialized to 5 kilometers
fRmax = CalculateRmax(fHmax);
if ( fDebugInfo ) { // if debug display results for all fith methods
fBetaInit = 2*ATan(CLUCE/(fAngularSpeedLF*fRmax)) - fData.fMedAlpha;
fBeta = PI - fBetaInit;
while(fBeta > 2*PI) fBeta = fBeta - 2*PI;
CalculatefromBetaEASdir();
Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirLF = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
fBetaInit = 2*ATan(CLUCE/(fAngularSpeedMF*fRmax)) - fData.fMedAlpha;
fBeta = PI - fBetaInit;
while(fBeta> 2*PI) fBeta = fBeta - 2*PI;
CalculatefromBetaEASdir();
Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirMF = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
fBetaInit = 2*ATan(CLUCE/(fAngularSpeedHF*fRmax)) - fData.fMedAlpha;
fBeta = PI - fBetaInit;
while(fBeta > 2*PI) fBeta = fBeta - 2*PI;
CalculatefromBetaEASdir();
Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirHF = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
}
// first value of beta
fBetaInit = 2*ATan(CLUCE/(fAngularSpeed*fRmax)) - fData.fMedAlpha;
fBeta = PI - fBetaInit;
while(fBeta > 2*PI) fBeta = fBeta - 2*PI;
CalculatefromBetaEASdir();
// recalculate Hmax
fHmax = FindHmax( fTHETAreco );
fRmax = CalculateRmax( fHmax );
TVector3 nmax = fData.fMedDir;
fData.VectorRmax.SetMagThetaPhi(fRmax,nmax.Theta(),nmax.Phi());
// recalculate beta
fBetaInit = 2*atan(CLUCE/(fAngularSpeed*fRmax)) - fData.fMedAlpha;
fBeta = PI - fBetaInit;
while(fBeta > 2*PI) fBeta = fBeta - 2*PI;
CalculatefromBetaEASdir();
fHmax = FindHmax(fTHETAreco);
Msg(EsafMsg::Info) << "method AA1(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
fAA1done = kTRUE;
}
//______________________________________________________________________________
void TrackDirection2Module::AA2() {
//
// ANALYTIC APPROXIMATED 2 method.
// Approximation: Shower velocity on a plane perpendicular to the detector
// axis is constant
//
Msg(EsafMsg::Info) << "Using method AA2()" << MsgDispatch;
if ( !fAA1done ) AA1();
if ( fAA1nan ) return;
vector<Double_t> xprov,yprov,errxyprov;
vector<Double_t> errxprov,erryprov;
for(Int_t i=0;i<fNumPoints;i++){
TVector3 dummy = fData.fSpPos[i];
xprov.push_back(dummy.x()*(fHISS-fHmax)/dummy.z());
yprov.push_back(dummy.y()*(fHISS-fHmax)/dummy.z());
Double_t dx = (fHISS-fHmax) * DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
Double_t dy = (fHISS-fHmax) * DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
errxyprov.push_back(1./Sqrt((Double_t)fData.fHits[i]));
errxprov.push_back( dx / fData.fHits[i]);
erryprov.push_back( dy / fData.fHits[i]);
}
Double_t speedx,speedy;
LeastSquaresFit *fitvx = new LeastSquaresFit(fNumPoints, fData.fTime, xprov, errxprov);
speedx = fitvx->GetSlope()*0.01; // why?
delete fitvx;
LeastSquaresFit *fitvy = new LeastSquaresFit(fNumPoints, fData.fTime, yprov, erryprov);
speedy = fitvy->GetSlope()*0.01; // why?
delete fitvy;
Double_t speed = Sqrt(speedx*speedx + speedy*speedy);
Double_t gthetamax = fData.fMedAlpha;
Double_t betacontrol = fBeta;
Double_t gammaV = 2*ATan((speed/CLUCE) * Sin(gthetamax));
Double_t betaV;
if(betacontrol*DEG>(90-gthetamax*DEG) && betacontrol*DEG<90){
betaV = -gammaV + gthetamax;
} else if (betacontrol*DEG > (90-gthetamax*DEG) && betacontrol*DEG > 90 ) {
betaV = gthetamax + gammaV;
}
fBeta = betaV;
while(fBeta > 2*PI) fBeta=fBeta-2*PI;
CalculatefromBetaEASdir();
fHmax = FindHmax(fTHETAreco);
Msg(EsafMsg::Info) << "method AA2(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
}
//______________________________________________________________________________
void TrackDirection2Module::NE1() {
//
// NUMERICAL EXACT 1 method
// Chi-square minimization of the difference between arrival times of photons
// measured and teoretically computed
//
Msg(EsafMsg::Info) << "using method NE1()" << MsgDispatch;
if( !fAA1done ) AA1();
if ( fAA1nan ) return;
TMinuit *minuitNE1 = new TMinuit(3);
minuitNE1->SetObjectFit(&fData);
minuitNE1->SetPrintLevel(fMinuitOutputLevel);
TString namebeta="beta", nametmax="tmax", namehmax="hmax";
Int_t ierflg(0);
Double_t fcnout,edm,errdef;
Int_t nvpar,nparx,icstat;
Double_t value, errv, vlow,vup;
Int_t gmaxcall=1000;
Double_t deltabeta,deltatmax;
Double_t betaStart = PI - fBeta;
Double_t hmaxStart = fHmax;
Double_t tmaxStart = fData.fMedTime*100;
Double_t step[3] = {0.2,2.500,0.3}; //about 10 deg (beta), 2500ns (tmax), 0.3 km (hmax)
Double_t arglist[10];
minuitNE1->SetFCN(ChiSquareNE1);
minuitNE1->mnparm(0, "beta",betaStart , step[0], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER beta.");
minuitNE1->mnparm(1, "tmax",tmaxStart , step[1], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER tmax.");
minuitNE1->mnparm(2, "hmax",hmaxStart , step[2], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER hmax.");
if( fFixTmaxNumeric ) {
arglist[0] = 2;
minuitNE1->mnexcm("FIX",arglist,1,ierflg); if (ierflg) Printf(" .....UNABLE to fix parameter tmax");
}
arglist[0] = 3;
minuitNE1->mnexcm("FIX",arglist,1,ierflg); if (ierflg) Printf(" .....UNABLE to fix parameter tmax");
arglist[0] = gmaxcall;
minuitNE1->mnexcm("MIGRAD",arglist ,1 , ierflg);if (ierflg) Printf(" ......UNABLE to minimize with MIGRAD");
minuitNE1->mnstat(fcnout,edm,errdef,nvpar,nparx,icstat);
minuitNE1->mnpout(0,namebeta, value,errv,vlow,vup,ierflg);
fBeta = PI - value;
deltabeta = errv;
minuitNE1->mnpout(2,nametmax, value,errv,vlow,vup,ierflg);
fTmaxFit = value;
deltatmax = errv;
CalculatefromBetaEASdir();
fHmax=FindHmax(fTHETAreco);
Msg(EsafMsg::Info) << "method NE1(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
}
//______________________________________________________________________________
void ChiSquareNE1(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag) {
//
// Function FCN called by minuit for the algorithm NE1
//
chisqnorm = 0;
ContainerData *fDataMin = (ContainerData*) gMinuit->GetObjectFit();
Int_t numhits(0);
Double_t chisq(0);
Double_t gerrort = fDataMin->fGtuLength;
Double_t Rm = (fDataMin->VectorRmax).Mag();
Double_t Am = fDataMin->fMedAlpha;
for(Int_t i=0;i<fDataMin->fNumPoints;i++){
Double_t TimeExpected = ( par[1] - Rm * ( ( Sin(Am-fDataMin->fAlpha[i]) + Sin(fDataMin->fAlpha[i]+par[0]) - Sin(Am+par[0]) )/ Sin(fDataMin->fAlpha[i]+par[0]) )/CLUCE);
chisq += Power((fDataMin->fTime[i]*100-TimeExpected)/gerrort,2)*fDataMin->fHits[i];
numhits += fDataMin->fHits[i];
}
chisqnorm = (Double_t)chisq/numhits;
}
//______________________________________________________________________________
void TrackDirection2Module::NE2() {
//
// NUMERICAL EXACT 2 method
// Chi-square minimization of angle between the versors of pixels in FOV
// and the corresponding vectors from points of the track and the detector
//
Msg(EsafMsg::Info) << "using method NE2()" << MsgDispatch;
if( !fAA1done ) AA1();
if ( fAA1nan ) return;
TMinuit *minuitNE2 = new TMinuit(3);
minuitNE2->SetObjectFit(&fData);
minuitNE2->SetPrintLevel(fMinuitOutputLevel);
TString nametheta="theta",namephi="phi",nametmax="tmax";
Int_t ierflg(0);
Double_t fcnout,edm,errdef;
Int_t nvpar,nparx,icstat;
Double_t value, errv, vlow,vup;
Int_t gmaxcall = 1000;
Double_t deltatheta,deltaphi,deltatmax;
Double_t thetastart = fTHETAreco;
Double_t phistart = fPHIreco;
Double_t tmaxStart = fData.fMedTime*100; //microseconds
Double_t step[3] = {0.2,0.2,2.500};
Double_t arglist[10];
minuitNE2->SetFCN(ChiSquareNE2);
minuitNE2->mnparm(0, "theta",thetastart , step[0], 0,0,ierflg); if (ierflg) Printf("......UNABLE TO DEFINE PARAMETER theta.");
minuitNE2->mnparm(1, "phi",phistart , step[1], 0,0,ierflg); if (ierflg) Printf(".....UNABLE TO DEFINE PARAMETER phi.");
minuitNE2->mnparm(2, "tmax",tmaxStart , step[2], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER tmax.");
if( fFixTmaxNumeric ) {
arglist[0] = 3;
minuitNE2->mnexcm("FIX",arglist,1,ierflg);if (ierflg) Printf(" .....UNABLE to fix parameter tmax");
}
arglist[0] = gmaxcall;
minuitNE2->mnexcm("MIGRAD",arglist ,1 , ierflg);if (ierflg) Printf(" .......UNABLE to minimize with migrad");
minuitNE2->mnstat(fcnout,edm,errdef,nvpar,nparx,icstat);
minuitNE2->mnpout(0,nametheta, value,errv,vlow,vup,ierflg);
fTHETAreco = value;
deltatheta = errv;
minuitNE2->mnpout(1,namephi, value,errv,vlow,vup,ierflg);
fPHIreco = value;
deltaphi = errv;
minuitNE2->mnpout(2,nametmax, value,errv,vlow,vup,ierflg);
fTmaxFit = value;
deltatmax = errv;
fEASDir.SetMagThetaPhi(1,fTHETAreco,fPHIreco);
fTHETAreco = fEASDir.Theta();
fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*PI) );
fDeltaTheta = fTHETAreco-fTrueTheta;
fDeltaPhi = fPHIreco-fTruePhi;
fDeltaEASDir = fEASDir.Angle(fTrueDir);
fHmax = FindHmax(fTHETAreco);
Msg(EsafMsg::Info) << "method NE2(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
}
//______________________________________________________________________
void ChiSquareNE2(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag) {
//
// Function FCN called by minuit for the algorithm NE2
//
ContainerData *fDataMin = (ContainerData*) gMinuit->GetObjectFit();
Int_t numhits(0);
Double_t PidotRi(0), modRi(0), Li(0), betaprimo(0);
TVector3 dummyRmax = fDataMin->VectorRmax;
Double_t Xm = dummyRmax.x();
Double_t Ym = dummyRmax.y();
Double_t Zm = dummyRmax.z();
Double_t Rm = dummyRmax.Mag();
Double_t gerrorpix = 0.0018; //rad, about 0.1 deg
Double_t sctmp = Sin(par[0])*Cos(par[1]-PI);
Double_t sstmp = Sin(par[0])*Sin(par[1]-PI);
Double_t ctmp = -Cos(par[0]);
Double_t thetaprimotmp = ACos((-Xm*sctmp-Ym*sstmp-Zm*ctmp)/Rm);
chisqnorm=0;
for(Int_t i=0; i<fDataMin->fNumPoints; i++){
betaprimo = 2*ATan(1/((1/Tan(thetaprimotmp/2))+CLUCE*(fDataMin->fTime[i]*100-par[2])/(Rm*Sin(thetaprimotmp))));
Li = Rm*(Cos(thetaprimotmp) - Sin(thetaprimotmp) / Tan(betaprimo));
PidotRi = fDataMin->fSpPos[i].x()*(Xm+Li*sctmp) + fDataMin->fSpPos[i].y()*(Ym+Li*sstmp) + fDataMin->fSpPos[i].z()*(Zm+Li*ctmp);
modRi = Power(Xm+Li*sctmp,2) + Power(Ym+Li*sstmp,2) + Power(Zm+Li*ctmp,2);
chisqnorm += Power((ACos(PidotRi/Sqrt(modRi)))/gerrorpix,2)*fDataMin->fHits[i];
numhits += fDataMin->fHits[i];
}
chisqnorm = (Double_t)chisqnorm/numhits;
}
//______________________________________________________________________________
void TrackDirection2Module::AE1() {
//
// ANALYTICAL EXACT 1 method
// Fit using exact relations between pixel directions in FOV and photons
// arrival times. This method doesn't require the knowledge of the TDP.
//
Msg(EsafMsg::Info) << "Using method AE1()" << MsgDispatch;
if( !fAA1done ) AA1();
if ( fAA1nan ) return;
vector<Double_t> alpha;
vector<Double_t> Ri;
vector<Double_t> dLi;
vector<Double_t> timestart;
vector<Double_t> nh;
vector<Double_t> x;
vector<Double_t> y;
vector<Double_t> z;
vector<Double_t> ex;
vector<Double_t> ey;
vector<Double_t> ez;
Double_t gnthetamax = fData.fMedDir.Theta();
Double_t gnphimax = fData.fMedDir.Phi();
Double_t timemax = fData.fMedTime*100; // microseconds
Int_t pointsAE1(0),hitsAE1(0);
Double_t Ritmp(0), dLitmp(0);
for ( Int_t i=0 ; i<fNumPoints; i++ ){
TVector3 dummy = fData.fSpPos[i];
Double_t dummyTime = fData.fTime[i]*100;//mus
Double_t cost = Cos(dummy.Theta())*Cos(gnthetamax) + Sin(dummy.Theta())*Sin(gnthetamax)*Cos(dummy.Phi()-gnphimax);
if( Abs(dummyTime-timemax)/fGtuLength>1 && Abs(dummyTime-timemax)/fGtuLength>fData.fApparentTimeLenght/10 && fData.fHits[i]>2 ){
if( dummyTime > timemax){
Ritmp = (2*CLUCE*(dummyTime-timemax)*fRmax + Power(CLUCE*(dummyTime-timemax),2))/(2*(fRmax*(1-cost)
+ CLUCE*(dummyTime-timemax)));
dLitmp = -Sqrt( Ritmp*Ritmp + fRmax*fRmax - 2*Ritmp*fRmax*cost );
}
if(dummyTime < timemax){
Ritmp = (2*CLUCE*(-dummyTime+timemax)*fRmax - Power(CLUCE*(dummyTime-timemax),2))/(2*(fRmax*(-1+cost)
+ CLUCE*(-dummyTime+timemax)));
dLitmp = Sqrt( Ritmp*Ritmp + fRmax*fRmax - 2*Ritmp*fRmax*cost );
}
alpha.push_back(fData.fAlpha[i]);
Ri.push_back(Ritmp);
dLi.push_back(dLitmp);
timestart.push_back(-dLitmp/CLUCE);
nh.push_back(fData.fHits[i]);
x.push_back(Ritmp*Sin(dummy.Theta())*Cos(dummy.Phi()));
ex.push_back(0.1/fData.fHits[i]);
y.push_back(Ritmp*Sin(dummy.Theta())*Sin(dummy.Phi()));
ey.push_back(0.1/fData.fHits[i]);
z.push_back(Ritmp*Cos(dummy.Theta()));
ez.push_back(0.1/fData.fHits[i]);
hitsAE1 += fData.fHits[i];
pointsAE1++;
}
}
alpha.push_back(fData.fMedAlpha);
dLi.push_back(0);
timestart.push_back(0);
Ri.push_back(fRmax);
x.push_back(fData.VectorRmax.x());
y.push_back(fData.VectorRmax.y());
z.push_back(fData.VectorRmax.z());
Double_t vx,vy,vz;
LeastSquaresFit *fitvx = new LeastSquaresFit(pointsAE1, timestart, x, ex);
vx = fitvx->GetSlope();
delete fitvx;
LeastSquaresFit *fitvz = new LeastSquaresFit(pointsAE1, timestart, z, ez);
vz = fitvz->GetSlope();
delete fitvz;
LeastSquaresFit *fitvy = new LeastSquaresFit(pointsAE1, timestart, y, ey);
vy = fitvy->GetSlope();
delete fitvy;
if ( IsNaN(vx) ||IsNaN(vy) || IsNaN(vz) ) {
Msg(EsafMsg::Warning) << "AE1() fit returns a NaN, skipping" << MsgDispatch;
return;
}
Double_t m1 = vx/vz;
Double_t m2 = vy/vz;
Double_t tT = Sqrt(m1*m1+m2*m2);
Double_t tP = m2/m1;
fTHETAreco = ATan(tT);
fPHIreco = ATan(tP);
if( tP>0 ){
fPHIreco = ( m2>0 ? atan(tP) : PI+atan(tP) );
} else {
fPHIreco = ( m2>0 ? PI+atan(tP) : (2*PI)+atan(tP) );
}
fEASDir.SetMagThetaPhi(1,fTHETAreco,fPHIreco);
fTHETAreco = fEASDir.Theta();
fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*PI) );
fDeltaTheta = fTHETAreco - fTrueTheta;
fDeltaPhi = fPHIreco - fTruePhi;
fDeltaEASDir = fEASDir.Angle(fTrueDir);
fHmax = FindHmax(fTHETAreco);
Msg(EsafMsg::Info) << "method AE1(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
}
//______________________________________________________________________________
void TrackDirection2Module::all() {
//
// Execute all methods for reconstructing the shower direction.
//
Msg(EsafMsg::Info) << "Using all methods .. " << MsgDispatch;
AA1();
fTHETArecoAA1 = fTHETAreco;
fPHIrecoAA1 = fPHIreco;
fDeltaEASDirAA1 = fDeltaEASDir;
AA2();
fTHETArecoAA2 = fTHETAreco;
fPHIrecoAA2 = fPHIreco;
fDeltaEASDirAA2 = fDeltaEASDir;
NE1();
fTHETArecoNE1 = fTHETAreco;
fPHIrecoNE1 = fPHIreco;
fDeltaEASDirNE1 = fDeltaEASDir;
NE2();
fTHETArecoNE2 = fTHETAreco;
fPHIrecoNE2 = fPHIreco;
fDeltaEASDirNE2 = fDeltaEASDir;
AE1();
fTHETArecoAE1 = fTHETAreco;
fPHIrecoAE1 = fPHIreco;
fDeltaEASDirAE1 = fDeltaEASDir;
//if all is chosen are saved only AA1 results
fTHETAreco = fTHETArecoAA1;
fPHIreco = fPHIrecoAA1;
}
//______________________________________________________________________________
void TrackDirection2Module::CalculatefromBetaEASdir() {
//
// Calculate vector of EAS direction using angle fBeta
// and the equation of TrackDirectionPlane
//
fEASDir = Cos(fBeta) * fW + Sin(fBeta) * fU;
fTHETAreco = fEASDir.Theta();
fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*PI) );
fDeltaTheta = fTHETAreco - fTrueTheta;
fDeltaPhi = fPHIreco - fTruePhi;
fDeltaEASDir = fEASDir.Angle(fTrueDir);
return;
}
//______________________________________________________________________________
Double_t TrackDirection2Module::FindHmax( Double_t theta ) {
//
// Find Hmax with the Linsley parametrization of the atmosphere
// depending on the zenith angle of the shower and the value of Xmax
// (fixed value is only a first approximation).
//
Double_t Xmax = 831; //g/cm^2
Double_t thetaloc = theta; //it should be zenith angle in local reference frame, to improve.
Double_t Xv = Xmax*Cos(thetaloc);
Double_t X0 = 1036.1; // g/cm^2
Double_t X4 = 631.1; // g/cm^2
Double_t X10 = 271.1; // g/cm^2
Double_t al(0), bl(0), cl(0);
Double_t val(0);
if( Xv<X0 && Xv>X4 ){
cl = 9.9418638; // km
bl = 1222.6562; // g/cm2
al = -186.5562; // g/cm2
} else if( Xv<X4 && Xv>X10 ) {
cl = 8.7815355; // km
bl = 1144.9069; // g/cm2
al = -94.9199; // g/cm2
} else if( Xv<X10 ) {
cl = 6.3614304; // km
bl = 1305.5948; // g/cm2
al = 0.61289; // g/cm2
} else {
Msg(EsafMsg::Warning) <<"...something was wrong in TrackDirection2Module::FindHmax"<< MsgDispatch;
return 0;
}
val = -cl * Log((Xmax * Cos(thetaloc)-al)/bl);
return val;
}
//______________________________________________________________________________
Double_t TrackDirection2Module::CalculateRmax( Double_t hmax ) {
//
// Method to geometrically find Rmax using Hmax and the versor pointing
// to the maximum of the shower
//
Double_t thmax = fData.fMedDir.Theta();
Double_t rmax = (RT+fHISS)*TMath::Cos(thmax) - TMath::Sqrt(TMath::Power(RT+hmax,2) - TMath::Power((RT+fHISS)*TMath::Sin(thmax),2));
return rmax;
}
//______________________________________________________________________________
Double_t TrackDirection2Module::DeltaX( Double_t theta, Double_t phi, Double_t errtheta, Double_t errphi ) {
//
// Calculate error on X projection on the focal surface of a given point
// on the unitary sphere
//
return Sqrt( Power((Cos(phi)*Cos(theta)+Sin(theta))*errtheta/(Cos(theta)*Cos(theta)),2.) +
Power(Sin(theta)*Sin(phi)*errphi/Cos(theta),2.));
}
//______________________________________________________________________________
Double_t TrackDirection2Module::DeltaY( Double_t theta, Double_t phi, Double_t errtheta, Double_t errphi ) {
//
// Calculate error on Y projection on the focal surface of a given point
// on the unitary sphere
//
return Sqrt( Power((Sin(phi)*Cos(theta)+Sin(theta))*errtheta/(Cos(theta)*Cos(theta)),2.) +
Power(Sin(theta)*Cos(phi)*errphi/Cos(theta),2.));
}
//___________________________________________________
void ContainerData::Clear() {
//
// Clear of the container data class
//
fPixXYZ.clear();
fPos.clear();
fErr.clear();
fSpPos.clear();
fSpErr.clear();
fSpPosSel.clear();
fTime.clear();
fTimeSel.clear();
fAlpha.clear();
fSigmaPhi.clear();
fSigmaTheta.clear();
fSigmaPhiSel.clear();
fSigmaThetaSel.clear();
fHits.clear();
fHitsSel.clear();
fNumPoints = 0;
fNumHits = 0;
fCentroid.SetXYZ(0,0,0);
fMedDir.SetXYZ(0,0,0);
VectorRmax.SetXYZ(0,0,0);
fMedTime = 0;
fMedAlpha = 0;
fApparentTimeLenght = 0;
fGtuLength = 0;
}