// $Id: TrackDirection2Module.cc,v 1.6 2005/04/14 12:38:58 moreggia Exp $
// Author: elena Nov, 19 2004
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: TrackDirection2Module *
* Package: <packagename> *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// TrackDirection2Module
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "TrackDirection2Module.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "LeastSquaresFit.hh"
#include "MedianFit.hh"
#include "RecoRootEvent.hh"
#include "EConst.hh"
#include "Config.hh"
#include <TGraph.h>
#include <TCanvas.h>
#include <TMinuit.h>
#include <TVector3.h>
Double_t DEG=TMath::RadToDeg();
Double_t PI=TMath::Pi();
Double_t RT=EConst::EarthRadius()/km;
Double_t CLUCE=EConst::C()*1e-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() {
Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
ConfigFileParser *pConfigRecoTrackDirection2Module = Config::Get()->GetCF("Reco","TrackDirection2Module");
fStat1 = pConfigRecoTrackDirection2Module->GetNum("TrackDirection2Module.fStat1");
fStat2 = pConfigRecoTrackDirection2Module->GetNum("TrackDirection2Module.fStat2");
fDoGraphUseShape = (Int_t)pConfigRecoTrackDirection2Module->GetNum("TrackDirection2Module.DoGraphUseShape");
fNumPointsMin = pConfigRecoTrackDirection2Module->GetNum("TrackDirection2Module.fNumPointsMin");
fNumHitsMinimum =(Int_t) pConfigRecoTrackDirection2Module->GetNum("TrackDirection2Module.fNumHitsMinimum");
fErrAngle = pConfigRecoTrackDirection2Module->GetNum("TrackDirection2Module.fErrAngle")/DEG;
if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fMethod")=="AA1") fMethodIdentifier = 1;
else if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fMethod")=="NE1") fMethodIdentifier = 2;
else if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fMethod")=="AA2") fMethodIdentifier = 3;
else if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fMethod")=="NE2") fMethodIdentifier = 4;
else if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fMethod")=="AE1") fMethodIdentifier = 5;
else if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fMethod")=="all") fMethodIdentifier = 6;
else {
Msg(EsafMsg::Panic) << "Wrong config value in TrackDirection2Module.cfg" << MsgDispatch;
throw runtime_error("Bad TrackDirection2Module.fMethod");
}
if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fUseShape")=="no") {
fDoShapeSelection = 0;
} else if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fUseShape")=="yes") {
fDoShapeSelection = 1;
fMulti1=(Int_t)pConfigRecoTrackDirection2Module->GetNum("TrackDirection2Module.fMulti1");
fMulti2=(Int_t)pConfigRecoTrackDirection2Module->GetNum("TrackDirection2Module.fMulti2");
} else {
Msg(EsafMsg::Panic) << " Wrong config value in TrackDirection2Module.cfg." << MsgDispatch;
throw runtime_error("Bad TrackDirection2Module.fUseShape");
}
if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fUseShapeInModules")=="yes") {
fUseShapeSelectioninModules = 1;
} else if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fUseShapeInModules")=="no") {
fUseShapeSelectioninModules = 0;
} else {
Msg(EsafMsg::Panic) << " Wrong config value in TrackDirection2Module.cfg." << MsgDispatch;
throw runtime_error("Bad TrackDirection2Module.fUseShapeInModules");
}
if(pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fFixTmaxNumeric")=="yes") {
fFixTmaxNumeric = 1;
} else if (pConfigRecoTrackDirection2Module->GetStr("TrackDirection2Module.fFixTmaxNumeric")=="no") {
fFixTmaxNumeric = 0;
} else {
Msg(EsafMsg::Panic) << " Wrong config value in TrackDirection2Module.cfg." << MsgDispatch;
throw runtime_error("Bad TrackDirection2Module.fFixTmaxNumeric");
}
ConfigFileParser *pConfigGeneralEuso = Config::Get()->GetCF("General","Euso");
HISS = pConfigGeneralEuso->GetNum("Euso.fAltitude");
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::PreProcess() {
fData.Clear();
fAA1done = 0;
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;
if(fDoGraphUseShape==1){
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");
if ( gcm==NULL ) {
Msg(EsafMsg::Warning) << " not found GTUClusteringModule." << MsgDispatch;
Int_t fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
Msg(EsafMsg::Info) <<"fTotalNumPoints="<<fTotalNumPoints<<MsgDispatch;
if ( fTotalNumPoints <= 0 ) {
Msg(EsafMsg::Info) << "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->GetObj("CluPixels") == NULL ) {
Msg(EsafMsg::Warning) << "No data in event" << MsgDispatch;
return kFALSE;
}
fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
Msg(EsafMsg::Info) << "pattern recognition data found " << MsgDispatch;
if (fPointsId.size() == 0) {
fNumPoints = 0;
Msg(EsafMsg::Warning) << "No pixels selected: terminated." << MsgDispatch;
return kFALSE;
} else {
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);
fData.fSigmaPhi.push_back(pix->GetPhiSigma());
fData.fSigmaTheta.push_back(pix->GetThetaSigma());
fData.fTime.push_back((Double_t)2.500*pix->GetGtu());
fData.fHits.push_back(pix->GetCounts());
fNumHits += pix->GetCounts();
}
fData.fNumPoints = fNumPoints;
fData.fNumHits = fNumHits;
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) << "FROM TRUTH : fTrueTheta = " << fTrueTheta*DEG << " deg,tfTruePhi = " << fTruePhi*DEG << " deg,tfTrueEnergy = " << fEv->GetHeader().GetTrueEnergy() << " MeV." << MsgDispatch;
TVector3 fTrueMaxOldsdr = fEv->GetHeader().GetTrueShowerMaxPos();
Double_t xmax=fTrueMaxOldsdr.x()/km;
Double_t ymax=fTrueMaxOldsdr.y()/km;
Double_t zmax=HISS-(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;
Int_t ShapeSelectionDone(0);
if ( fDoShapeSelection == 0 ) {
FindPlane();
ShapeSelectionDone=0;
} else if ( fDoShapeSelection == 1 ) {
UseShapeandFindPlane();
if( fNumPointsSel < fNumPointsMin ){
Msg(EsafMsg::Info) << "Impossible to use selection by shape: only " << fNumPointsSel << " pixels selected!" <<MsgDispatch;
ShapeSelectionDone = 0;
FindPlane();
}else{ ShapeSelectionDone = 1; }
}
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( fUseShapeSelectioninModules == 1 && ShapeSelectionDone == 1 ){
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]);
}
}
if ( fMethodIdentifier == 1 ) AA1();
else if ( fMethodIdentifier == 2 ) NE1();
else if ( fMethodIdentifier == 3 ) AA2();
else if ( fMethodIdentifier == 4 ) NE2();
else if ( fMethodIdentifier == 5 ) AE1();
else if ( fMethodIdentifier == 6 ) all();
else {
Msg(EsafMsg::Panic) << "Wrong config value in TrackDirection2Module.cfg" << MsgDispatch;
throw runtime_error("Bad TrackDirection2Module.fMethod");
}
fData.Clear();
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::PostProcess() {
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() {
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;
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::DoStat(vector<Double_t> err, string namemethod){
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++;
}
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;
break;
}
}
string titleh = namemethod.c_str();
titleh += Form(" Resolution(68%)=%g deg", resolution);
hR->SetTitle(titleh.c_str());
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 kTRUE;
}
void TrackDirection2Module::UserMemoryClean() {
}
//_____________________________________________________________________________
Bool_t TrackDirection2Module::SaveRootData(RecoRootEvent *fRecoRootEvent) {
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::FindPlane(){
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;
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
Double_t p = dummy.Phi();
Double_t t = dummy.Theta();
Double_t dx = TMath::Sqrt( TMath::Power((cos(p)*cos(t)+sin(t))*fData.fSigmaTheta[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*sin(p)*fData.fSigmaPhi[i]/cos(t),2.));
Double_t dy = TMath::Sqrt( TMath::Power((sin(p)*cos(t)+sin(t))*fData.fSigmaTheta[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*cos(p)*fData.fSigmaPhi[i]/cos(t),2.));
xpro.push_back(dummy.X()/dummy.Z());
ypro.push_back(dummy.Y()/dummy.Z());
tpro.push_back(fData.fTime[i]);
errxpro.push_back(dx/fData.fHits[i]);
errypro.push_back(dy/fData.fHits[i]);
}
Msg(EsafMsg::Info) <<"using LeastSquaresFit in FindPlane()"<< MsgDispatch;
LeastSquaresFit *xtFit = new LeastSquaresFit(fNumPoints, tpro, xpro, errxpro);
Double_t vX = xtFit->GetSlope();
Double_t wX = xtFit->GetOffset();
delete xtFit;
LeastSquaresFit *ytFit = new LeastSquaresFit(fNumPoints, tpro, ypro, errypro);
Double_t vY = ytFit->GetSlope();
Double_t wY = ytFit->GetOffset();
delete ytFit;
fNorm.SetXYZ(vY,-vX,wY*vX-wX*vY);
fNorm.SetMag(1.);
if ( fNorm.Z() > 0 ) fNorm *= -1;
xpro.clear();
ypro.clear();
tpro.clear();
errxpro.clear();
errypro.clear();
}
//______________________________________________________________________________
void TrackDirection2Module::UseShapeandFindPlane(){
vector<Double_t> xpro;
vector<Double_t> ypro;
vector<Double_t> errxpro;
vector<Double_t> errypro;
vector<Double_t> tpro;
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
Double_t p = dummy.Phi();
Double_t t = dummy.Theta();
Double_t dx = TMath::Sqrt( TMath::Power((cos(p)*cos(t)+sin(t))*fData.fSigmaTheta[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*sin(p)*fData.fSigmaPhi[i]/cos(t),2.));
Double_t dy = TMath::Sqrt( TMath::Power((sin(p)*cos(t)+sin(t))*fData.fSigmaTheta[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*cos(p)*fData.fSigmaPhi[i]/cos(t),2.));
xpro.push_back(dummy.X()/dummy.Z());
ypro.push_back(dummy.Y()/dummy.Z());
tpro.push_back(fData.fTime[i]);
errxpro.push_back(dx/fData.fHits[i]);
errypro.push_back(dy/fData.fHits[i]);
}
MedianFit *xMf1 = new MedianFit(fNumPoints, tpro, xpro);
Double_t sX = xMf1->GetSlope();
Double_t oX = xMf1->GetOffset();
Double_t aX = xMf1->GetAbsoluteDeviation();
delete xMf1;
MedianFit *yMf1 = new MedianFit(fNumPoints, tpro, ypro);
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==1){
Double_t xprog[fNumPoints],yprog[fNumPoints],tprog[fNumPoints],xretta[fNumPoints],yretta[fNumPoints];
Double_t xrettaup[fNumPoints],yrettaup[fNumPoints];
Double_t xrettadown[fNumPoints],yrettadown[fNumPoints];
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
xretta[i]=sX*fData.fTime[i]+oX;
xrettaup[i]=xretta[i]+aX;
xrettadown[i]=xretta[i]-aX;
yretta[i]=sY*fData.fTime[i]+oY;
yrettaup[i]=yretta[i]+aY;
yrettadown[i]=yretta[i]-aY;
xprog[i]=dummy.X()/dummy.Z();
yprog[i]=dummy.Y()/dummy.Z();
tprog[i]=fData.fTime[i];
}
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,xretta);
gxproretta->SetMarkerColor(4);gxproretta->SetLineColor(4);
TGraph *gyproretta = new TGraph(fNumPoints,tprog,yretta);
gyproretta->SetMarkerColor(4);gyproretta->SetLineColor(4);
TGraph *gxprorettaup = new TGraph(fNumPoints,tprog,xrettaup);
gxprorettaup->SetMarkerColor(2);gxprorettaup->SetLineColor(2);
TGraph *gyprorettaup = new TGraph(fNumPoints,tprog,yrettaup);
gyprorettaup->SetMarkerColor(2);gyprorettaup->SetLineColor(2);
TGraph *gxprorettadown = new TGraph(fNumPoints,tprog,xrettadown);
gxprorettadown->SetMarkerColor(2);gxprorettadown->SetLineColor(2);
TGraph *gyprorettadown = new TGraph(fNumPoints,tprog,yrettadown);
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();
}
//////////////////////////////////////////////////////////////////////////////
vector<Double_t> xprosel;
vector<Double_t> yprosel;
vector<Double_t> errxprosel;
vector<Double_t> erryprosel;
vector<Double_t> tprosel;
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];
if(TMath::Abs(xprotmp-sX*tprotmp-oX)<fMulti1*aX || TMath::Abs(yprotmp-sY*tprotmp-oY)<fMulti1*aY){
Double_t p = dummy.Phi();
Double_t t = dummy.Theta();
Double_t dx = TMath::Sqrt( TMath::Power((cos(p)*cos(t)+sin(t))*fData.fSigmaTheta[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*sin(p)*fData.fSigmaPhi[i]/cos(t),2.));
Double_t dy = TMath::Sqrt( TMath::Power((sin(p)*cos(t)+sin(t))*fData.fSigmaTheta[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*cos(p)*fData.fSigmaPhi[i]/cos(t),2.));
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];
}
}
Msg(EsafMsg::Info) << "UseShapeandFindPlane(), first step: selected " << fNumPointsSeltmp << " from " << fNumPoints << " pixels" << MsgDispatch;
if ( fNumPointsSeltmp < fNumPointsMin ){
fNumPointsSel = fNumPointsSeltmp;
return;
}
MedianFit *xMf2 = new MedianFit(fNumPointsSeltmp, tprosel, xprosel);
Double_t sX2 = xMf2->GetSlope();
Double_t oX2 = xMf2->GetOffset();
Double_t aX2 = xMf2->GetAbsoluteDeviation();
delete xMf2;
MedianFit *yMf2 = new MedianFit(fNumPointsSeltmp, tprosel, yprosel);
Double_t sY2 = yMf2->GetSlope();
Double_t oY2 = yMf2->GetOffset();
Double_t aY2 = yMf2->GetAbsoluteDeviation();
delete yMf2;
vector<Double_t> xprosel2;
vector<Double_t> yprosel2;
vector<Double_t> errxprosel2;
vector<Double_t> erryprosel2;
vector<Double_t> tprosel2;
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];
if(TMath::Abs(xprotmp-sX2*tprotmp-oX2)<fMulti2*aX2 || TMath::Abs(yprotmp-sY2*tprotmp-oY2)<fMulti2*aY2){
Double_t p = dummy.Phi();
Double_t t = dummy.Theta();
Double_t dx = TMath::Sqrt( TMath::Power((cos(p)*cos(t)+sin(t))*fSigmaThetaSeltmp[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*sin(p)*fSigmaPhiSeltmp[i]/cos(t),2.));
Double_t dy = TMath::Sqrt( TMath::Power((sin(p)*cos(t)+sin(t))*fSigmaThetaSeltmp[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*cos(p)*fSigmaPhiSeltmp[i]/cos(t),2.));
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]);
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];
}
}
Msg(EsafMsg::Info) << "UseShapeandFindPlane(), second step: selected " << fNumPointsSel << " from " << fNumPointsSeltmp << " pixels" << MsgDispatch;
if ( fNumPointsSel < fNumPointsMin ) return;
/*LeastSquaresFit *xtFit = new LeastSquaresFit(fNumPointsSel, tprosel, xprosel, errxprosel);
Double_t vX = xtFit->GetSlope();
Double_t wX = xtFit->GetOffset();
delete xtFit;
LeastSquaresFit *ytFit = new LeastSquaresFit(fNumPointsSel, tprosel, yprosel, erryprosel);
DoMedianFit *xtFit = new MedianFit(fNumPointsSel, tprosel, xprosel);
Double_t vX = xtFit->GetSlope();
Double_t wX = xtFit->GetOffset();
delete xtFit;*/
MedianFit *xtFit = new MedianFit(fNumPointsSel, tprosel, xprosel);
Double_t vX = xtFit->GetSlope();
Double_t wX = xtFit->GetOffset();
Double_t a2X = xtFit->GetAbsoluteDeviation();
delete xtFit;
MedianFit *ytFit = new MedianFit(fNumPointsSel, tprosel, yprosel);
Double_t vY = ytFit->GetSlope();
Double_t wY = ytFit->GetOffset();
Double_t a2Y = ytFit->GetAbsoluteDeviation();
delete ytFit;
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==1){
Double_t xprog2[fNumPointsSel],yprog2[fNumPointsSel],tprog2[fNumPointsSel],xretta2[fNumPointsSel],yretta2[fNumPointsSel];
Double_t xretta2up[fNumPoints],yretta2up[fNumPoints];
Double_t xretta2down[fNumPoints],yretta2down[fNumPoints];
for( Int_t i=0; i<fNumPointsSel; i++ ) {
TVector3 dummy = fData.fSpPosSel[i];
xretta2[i]=vX*fData.fTimeSel[i]+wX;
yretta2[i]=vY*fData.fTimeSel[i]+wY;
xretta2up[i]=xretta2[i]+a2X;
xretta2down[i]=xretta2[i]-a2X;
yretta2up[i]=yretta2[i]+a2Y;
yretta2down[i]=yretta2[i]-a2Y;
xprog2[i]=dummy.X()/dummy.Z();
yprog2[i]=dummy.Y()/dummy.Z();
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,xretta2);
gxproretta2->SetMarkerColor(4);gxproretta2->SetLineColor(4);
TGraph *gyproretta2 = new TGraph(fNumPointsSel,tprog2,yretta2);
gyproretta2->SetMarkerColor(4);gyproretta2->SetLineColor(4);
TGraph *gxproretta2up = new TGraph(fNumPointsSel,tprog2,xretta2up);
gxproretta2up->SetMarkerColor(2);gxproretta2up->SetLineColor(2);
TGraph *gyproretta2up = new TGraph(fNumPointsSel,tprog2,yretta2up);
gyproretta2up->SetMarkerColor(2);gyproretta2up->SetLineColor(2);
TGraph *gxproretta2down = new TGraph(fNumPointsSel,tprog2,xretta2down);
gxproretta2down->SetMarkerColor(2);gxproretta2down->SetLineColor(2);
TGraph *gyproretta2down = new TGraph(fNumPointsSel,tprog2,yretta2down);
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();
}
//////////////////////////////////////////////////////////////////////////////
xpro.clear();
ypro.clear();
tpro.clear();
errxpro.clear();
errypro.clear();
xprosel.clear();
yprosel.clear();
tprosel.clear();
errxprosel.clear();
erryprosel.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() {
Msg(EsafMsg::Info) << "using method AA1()" << MsgDispatch;
fData.fMedDir.SetXYZ(0,0,0);
fData.fMedTime=0;
fData.fMedAlpha=0;
vector<Double_t> erralpha;
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]));
fData.fMedTime += fData.fTime[i]*fData.fHits[i];
fData.fMedDir += fData.fSpPos[i]*fData.fHits[i];
}
fData.fApparentTimeLenght = (Int_t)((fData.fTime[fNumPoints-1]-fData.fTime[0])/2.500);
fData.fMedTime = fData.fMedTime/fNumHits;
fData.fMedDir = (fData.fMedDir*(1/(Double_t) fNumHits)).Unit();
fData.fMedAlpha = acos(fData.fMedDir*fW);
//use the following only to debug (plot angles alphai vs time)
/*Double_t alphai[10000],timei[10000];
for( Int_t i=0; i<fNumPoints; i++ ) {
TVector3 dummy = fData.fSpPos[i];
alphai[i]=acos(dummy*fW)*DEG;
timei[i]=fData.fTime[i];
}
TGraph *gang = new TGraph(fNumPoints,timei,alphai);
gang->SetNameTitle("gang","alphai vs time");
gang->SetMarkerSize(.5);
gang->SetMarkerStyle(23);
TCanvas *cAA1=new TCanvas("cAA1","cAA1",1);cAA1->Divide(1,1);
cAA1->cd(1);gang->Draw("AP");
cAA1->Write();
*/
LeastSquaresFit *ASl = new LeastSquaresFit(fNumPoints, fData.fTime, fData.fAlpha, erralpha);
fAngularSpeed = ASl->GetSlope();
/*
MedianFit *ASm = new MedianFit(fNumPoints, fData.fTime, fData.fAlpha);
fAngularSpeed = ASm->GetSlope();
*/
fHmax=5;
fRmax=CalculateRmax(fHmax);
fBetaInit = 2*atan(CLUCE/(fAngularSpeed*fRmax)) - fData.fMedAlpha;
fBeta = PI - fBetaInit;
CalculatefromBetaEASdir();
fHmax=FindHmax(fTHETAreco);
fRmax=CalculateRmax(fHmax);
TVector3 nmax=fData.fMedDir;
fData.VectorRmax.SetMagThetaPhi(fRmax,nmax.Theta(),nmax.Phi());
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=1;
}
//______________________________________________________________________________
void TrackDirection2Module::AA2() {
Msg(EsafMsg::Info) << "using method AA2()" << MsgDispatch;
if(!fAA1done) AA1();
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()*(HISS-fHmax)/dummy.z());
yprov.push_back(dummy.y()*(HISS-fHmax)/dummy.z());
Double_t p = dummy.Phi();
Double_t t = dummy.Theta();
Double_t dx = (HISS-fHmax) * TMath::Sqrt( TMath::Power((cos(p)*cos(t)+sin(t))*fData.fSigmaTheta[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*sin(p)*fData.fSigmaPhi[i]/cos(t),2.));
Double_t dy = (HISS-fHmax) * TMath::Sqrt( TMath::Power((sin(p)*cos(t)+sin(t))*fData.fSigmaTheta[i]/(cos(t)*cos(t)),2.)
+TMath::Power(sin(t)*cos(p)*fData.fSigmaPhi[i]/cos(t),2.));
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();
delete fitvx;
LeastSquaresFit *fitvy = new LeastSquaresFit(fNumPoints, fData.fTime, yprov, erryprov);
speedy = fitvy->GetSlope();
delete fitvy;
Double_t speed=sqrt(speedx*speedx+speedy*speedy);
Double_t gtetamax = fData.fMedAlpha;
Double_t betacontrol=fBeta;
Double_t gammaV=2*atan((speed/CLUCE)*sin(gtetamax));
Double_t betaV;
if(betacontrol*DEG>(90-gtetamax*DEG) && betacontrol*DEG<90){
betaV=-gammaV+gtetamax;
}else{
if(betacontrol*DEG>(90-gtetamax*DEG) && betacontrol*DEG>90){
betaV=gtetamax+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() {
Msg(EsafMsg::Info) << "using method NE1()" << MsgDispatch;
if(!fAA1done) AA1();
TMinuit *minuitNE1 = new TMinuit(3);
minuitNE1->SetObjectFit(&fData);
minuitNE1->SetPrintLevel(-1);
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;
Double_t step[3]={0.2,2.500,0.3};//about 10 deg (beta), 2500ns (tmax), 0.3 km (hmax)
Double_t arglist[10];
//Msg(EsafMsg::Info) << "NE1: initializing minuit fit with beta = "<<betaStart*DEG<<" deg" << MsgDispatch;
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==1 ) {
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;
//Msg(EsafMsg::Info) << "fBeta="<<fBeta*DEG<<" deg con deltabeta="<<deltabeta*DEG<<" deg da NE1"<< MsgDispatch;
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;
}
//**********Function FCN called by minuit for the algorithm NE1
void ChiSquareNE1(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag){
chisqnorm=0;
ContainerData *fDataMin = (ContainerData*) gMinuit->GetObjectFit();
Int_t numhits(0);
Double_t chisq(0);
Double_t gerrort=2.500;
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 += pow((fDataMin->fTime[i]-TimeExpected)/gerrort,2)*fDataMin->fHits[i];
numhits += fDataMin->fHits[i];
}
chisqnorm = (Double_t)chisq/numhits;
//cout<<"tmax =par[1]="<<par[1]<<" msthmax="<<par[2]<<" kmtbeta=par[0]="<<par[0]*DEG<<"tchisqnorm="<<chisqnorm<<endl;
}
//______________________________________________________________________________
void TrackDirection2Module::NE2() {
Msg(EsafMsg::Info) << "using method NE2()" << MsgDispatch;
if(!fAA1done) AA1();
TMinuit *minuitNE2 = new TMinuit(3);
minuitNE2->SetObjectFit(&fData);
minuitNE2->SetPrintLevel(-1);
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;
Double_t step[3]={0.2,0.2,2.500};
Double_t arglist[10];
minuitNE2->SetFCN(ChiSquareNE2);
//Msg(EsafMsg::Info) << "NE2: initializing minuit fit with theta = "<<thetastart*DEG<<" deg, phi = "<<phistart*DEG<<" deg" << MsgDispatch;
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==1 ) {
arglist[0] = 2;
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;
}
//***********Function FCN called by minuit for the algorithm NE2
void ChiSquareNE2(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag){
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]-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=pow(Xm+Li*sctmp,2)+pow(Ym+Li*sstmp,2)+pow(Zm+Li*ctmp,2);
chisqnorm += pow((acos(PidotRi/sqrt(modRi)))/gerrorpix,2)*fDataMin->fHits[i];
numhits += fDataMin->fHits[i];
}
chisqnorm = (Double_t)chisqnorm/numhits;
//cout <<"chisqnorm="<<chisqnorm<<"ttmax = par[2]="<<par[2]<<"ttheta=par[0]="<<par[0]*DEG<<"tphi=par[1]="<<par[1]*DEG<<endl;
}
//______________________________________________________________________________
void TrackDirection2Module::AE1() {
Msg(EsafMsg::Info) << "using method AE1()" << MsgDispatch;
if(!fAA1done) AA1();
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;
Int_t pointsNE1(0),hitsNE1(0);
Double_t Ritmp,dLitmp;
for ( Int_t i=0 ; i<fNumPoints; i++ ){
TVector3 dummy = fData.fSpPos[i];
Double_t dummyTime = fData.fTime[i];
Double_t cost=cos(dummy.Theta())*cos(gnthetamax)+sin(dummy.Theta())*sin(gnthetamax)*cos(dummy.Phi()-gnphimax);
if( TMath::Abs(dummyTime-timemax)/2.500 > 1 && TMath::Abs(dummyTime-timemax)/2.500 > fData.fApparentTimeLenght/10 && fData.fHits[i]>2 ){
if(dummyTime>timemax){
Ritmp = (2*CLUCE*(dummyTime-timemax)*fRmax + pow(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 - pow(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]);
hitsNE1 += fData.fHits[i];
pointsNE1++;
}
}
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(pointsNE1, timestart, x, ex);
vx = fitvx->GetSlope();
delete fitvx;
LeastSquaresFit *fitvz = new LeastSquaresFit(pointsNE1, timestart, z, ez);
vz = fitvz->GetSlope();
delete fitvz;
LeastSquaresFit *fitvy = new LeastSquaresFit(pointsNE1, timestart, y, ey);
vy = fitvy->GetSlope();
delete fitvy;
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() {
Msg(EsafMsg::Info) << ".. using all methods .. " << MsgDispatch;
if(!fAA1done) 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;
}
//______________________________________________________________________________
//method to calculate vector of EAS direction using angle fBeta and the equation of TrackDirectionPlane
void TrackDirection2Module::CalculatefromBetaEASdir() {
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;
}
//______________________________________________________________________________
//method to 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 TrackDirection2Module::FindHmax( Double_t theta ) {
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,bl,cl;
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*TMath::Log((Xmax*cos(thetaloc)-al)/bl);
//Msg(EsafMsg::Info) <<"FindHmax(): Xmax="<<Xmax<<" g/cm^2, Xv="<<Xv<<" g/cm^2, thetaloc="<<thetaloc*DEG <<" deg -----> Hmax = "<<val<<" km"<<MsgDispatch;
return val;
}
//______________________________________________________________________________
//method to geometrically find Rmax using Hmax and the versor pointing to the maximum of the shower
Double_t TrackDirection2Module::CalculateRmax( Double_t hmax ) {
Double_t thmax = fData.fMedDir.Theta();
Double_t rmax = (RT+HISS)*TMath::Cos(thmax) - TMath::Sqrt(TMath::Power(RT+hmax,2) - TMath::Power((RT+HISS)*TMath::Sin(thmax),2));
return rmax;
}
//___________________________________________________
void ContainerData::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;
}