Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TrackDirection2Module - source file

// $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;
	
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31