Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TrackDirection2Module - source file

// $Id: TrackDirection2Module.cc,v 1.16 2005/10/21 13:30:08 moreggia Exp $
// Author: elena   Nov, 19 2004

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: TrackDirection2Module                                                *
 *  Package: Fitting                                                         *
 *  Coordinator: <coordinator>                                               *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// TrackDirection2Module
//
// This module is devoted to the reconstruction of the shower direction.
// It implements some different algorithms (descibed in detail in ....)
// Essentially the module keeps the points found by pattern recogniton
// (clustering or Hough transform) and first find the plane that contains the 
// track and the detector (TDP).
// The reserarch of TDP can be done in the following ways:
//     - further selection of points with hough transform
//     - further selection of points with shape selection method
//     - no further selection
// In each case the TDP is founded by a fit (least squares, median or hough)
// of the x-t, y-t projections of points on the plane z=0.
//
// Then the shower direction is reconstructed by one of the following methods:
//     - analytical approximated 1 AA1()
//     - analytical approximated 2 AA2()
//     - numerical exact 1 NE1()
//     - numerical exact 2 NE2()
//     - analytical exact 1 AE1()
//
// The AA1() method is in each case used in order to initialize the fit for all
// other methods.
//
//   Config file parameters
//   ======================
//
//   fNumHitsMinimum : minimum number of hits per pixel in a GTU
//   fNumPointsMinimum : minimum number of points to proceed with reconstruction
//   fUseHough [bool] : use hough trasnform to find the TDP 
//   fUseHoughwithselection [bool] : use the points selected with Hough transform
//                                   in the reconstruction of direction
//   fUseShape [bool] : use shape selection method to find the TDP
//   fUseShapeInModules [bool] : use the points selected with shape selction
//                               in the reconstruction of direction
//   fDebugInfo [bool] : compute and display some debug informations
//   
//   fMulti1 : multiplier for shape selection method
//   fMulti2 : multiplier for shape selection method
//
//   fAA1FitMethod : fit method for AA1 algorithm
//   - Valid options : linear (least squares fit)
//                     median (median fit)
//                     hough  (hough fit)
//
//   fMethod : direction reconstruction method
//   - Valid options : AA1 || AA2 || NE1 || NE2 || AE1 (single algorithm)
//                     all                             (executes all algorithms)
//
//   fFixTmaxNumeric [bool] : fix shower maximum parameters in numerical fits
//   fErrAngle : angular error [deg]
//
//   fStat1 : angular error [deg] value 1 for statistics
//   fStat2 : angular error [deg] value 2 for statistics
//   - Statistics are made between 0 < err < fStat1 and fStat1 < err < fStat2
//
//   fDoGraphUseShape [bool] : save some debug graph in rootfile
//
//   fMinuitOutputLevel : set the MINUIT output display level
//

#include "EusoCluster.hh"
#include "TrackDirection2Module.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "LeastSquaresFit.hh"
#include "MedianFit.hh"
#include "HoughFit.hh"
#include "RecoRootEvent.hh"
#include "EConst.hh"
#include "Config.hh"
#include "ERunParameters.hh"

#include <TGraph.h>
#include <TCanvas.h>
#include <TMinuit.h>
#include <TVector3.h>
#include <TH2F.h>
#include <TLegend.h>

using namespace TMath;
using namespace sou;
using namespace EConst;

Double_t DEG = RadToDeg();
Double_t PI = Pi();
Double_t RT = EarthRadius()/km; //earth radius in km
Double_t CLUCE = Clight()*1.e-3;     //constant c in unit km/microsecond

//FCN function called by numeric method NE1
void ChiSquareNE1(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag);

//FCN function called by numeric method NE2
void ChiSquareNE2(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag);

ClassImp(TrackDirection2Module)
ClassImp(ContainerData)

//_____________________________________________________________________________
 TrackDirection2Module::TrackDirection2Module() : RecoModule("TrackDirection2") {
    //
    // ctor
    //
}

//_____________________________________________________________________________
 TrackDirection2Module::~TrackDirection2Module() {
    //
    // dtor
    //
}

//_____________________________________________________________________________
 Bool_t TrackDirection2Module::Init() {
    //
    // Initialization of variables
    //
	 
    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;

    useLFaa1 = kFALSE; useMFaa1 = kFALSE; useHFaa1 = kFALSE;
    useLFplane = kFALSE; useMFplane = kFALSE; useHFplane = kFALSE;
	
    fStat1 = Conf()->GetNum("TrackDirection2Module.fStat1");
    fStat2 = Conf()->GetNum("TrackDirection2Module.fStat2");
    fDoGraphUseShape = Conf()->GetBool("TrackDirection2Module.fDoGraphUseShape");
    fNumPointsMin = (Int_t)Conf()->GetNum("TrackDirection2Module.fNumPointsMin");
    fNumHitsMinimum = (Int_t)Conf()->GetNum("TrackDirection2Module.fNumHitsMinimum");
    fErrAngle = Conf()->GetNum("TrackDirection2Module.fErrAngle")/DEG;

    if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="linear") useLFplane = kTRUE;
    else if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="median") useMFplane = kTRUE;
    else if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="hough") useHFplane = kTRUE;
    else Msg(EsafMsg::Panic) << "Wrong config fTDPFitMethod value in TrackDirection2Module.cfg" << MsgDispatch; 
	
    if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="linear") useLFaa1 = kTRUE;
    else if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="median") useMFaa1 = kTRUE;
    else if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="hough") useHFaa1 = kTRUE;
    else Msg(EsafMsg::Panic) << "Wrong config fAA1FitMethod value in TrackDirection2Module.cfg" << MsgDispatch; 
	
    if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AA1")      fMethodIdentifier = 1;
    else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="NE1") fMethodIdentifier = 2;
    else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AA2") fMethodIdentifier = 3;
    else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="NE2") fMethodIdentifier = 4;
    else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AE1") fMethodIdentifier = 5;
    else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="all") fMethodIdentifier = 6;
    else Msg(EsafMsg::Panic) << "Wrong config fMethod value in TrackDirection2Module.cfg" << MsgDispatch;
	
    fDoHough = Conf()->GetBool("TrackDirection2Module.fUseHough");
    fOptionSelectionHough = Conf()->GetBool("TrackDirection2Module.fUseHoughwithselection");
    fDoShapeSelection = Conf()->GetBool("TrackDirection2Module.fUseShape");
    if ( fDoShapeSelection) {
	fMulti1 = Conf()->GetNum("TrackDirection2Module.fMulti1");
	fMulti2 = Conf()->GetNum("TrackDirection2Module.fMulti2");
    }
	
    fUseShapeSelectioninModules = Conf()->GetBool("TrackDirection2Module.fUseShapeInModules");
    fFixTmaxNumeric = Conf()->GetBool("TrackDirection2Module.fFixTmaxNumeric");
    fDebugInfo = Conf()->GetBool("TrackDirection2Module.fDebugInfo");
    fMinuitOutputLevel = (Int_t)Conf()->GetNum("TrackDirection2Module.fMinuitOutputLevel");
    if ( fMinuitOutputLevel < -1 || fMinuitOutputLevel > 3)
        Msg(EsafMsg::Panic) << "Wrong fMinuitOutputLevel value in TrackDirection2Module.cfg" << MsgDispatch;
    
    ConfigFileParser *pConfigGeneralEuso = Config::Get()->GetCF("General","Euso");
    fHISS = pConfigGeneralEuso->GetNum("Euso.fAltitude");

    // gtu length in microseconds
    ConfigFileParser *pConfigMacroCell = Config::Get()->GetCF("Electronics","MacroCell");
    fGtuLength = pConfigMacroCell->GetNum("MacroCell.fGtuTimeLength")/microsecond;
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t TrackDirection2Module::PreProcess() {
    //
    // Pre-process of the reco event
    //
     
    fData.Clear();
    fAA1done = kFALSE;
    fAA1nan = kFALSE;
    fNumPoints = 0;
    fNumHits = 0;
    fNumPointsSel = 0;
    fNumHitsSel = 0;
    fQuality = -1;

    fCentroid.SetXYZ(0,0,0);
    fNorm.SetXYZ(0,0,0);
    fNormsel.SetXYZ(0,0,0);
    fW.SetXYZ(0,0,0);
    fU.SetXYZ(0,0,0);
    fTrueDir.SetXYZ(0,0,0);
    fTrueNorm.SetXYZ(0,0,0);
    fTrueMax.SetXYZ(0,0,0);
    fEASDir.SetXYZ(0,0,0);

    fAngularSpeed = 0;
    fBeta = 0;
    fBetaInit = 0;
    fHmax = 0;
    fRmax = 0;
    fTrueTheta = 0;
    fTruePhi = 0;
    fTHETAloc = 0;
    fTHETAreco = 0;
    fPHIreco = 0;
    fTmaxFit = 0;
    fDeltaTheta = 0;
    fDeltaPhi = 0;
    fDeltaEASDir = 0;
    fDeltaTDP = 0;
 
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t TrackDirection2Module::Process(RecoEvent *ev) {
	
    fEv = ev;
    RecoEventHeader header = fEv->GetHeader();
    ERunParameters *runpars = header.GetRunPars();
    
    if(fDoGraphUseShape){
	gDirectory->cd("/");
    	string pathname = Form("TD2ProcessRecoEvent%d", fEv->GetHeader().GetNum());
    	TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
    	path->cd();
    }
    
	
    const RecoModuleData *gcm = fEv->GetModuleData("GTUClustering");
    const RecoModuleData *moduledataHough = fEv->GetModuleData("HoughTransform");
    if ( gcm == NULL && moduledataHough == NULL ) {

        Msg(EsafMsg::Warning) << "Not found GTUClusteringModule and not found HoughTransformModule" << MsgDispatch;
	    
        Int_t fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
        Msg(EsafMsg::Info) <<"fTotalNumPoints = "<<fTotalNumPoints<<MsgDispatch;
        if ( fTotalNumPoints <= 0 ) {       
            Msg(EsafMsg::Warning) << "No points in event! Exit from this event." << MsgDispatch;
            return kFALSE;
        }
        fPointsId.clear();
        fNumPoints = 0;
        for( Int_t i=0; i<fTotalNumPoints; i++ ) {
    	    if ( fEv->GetRecoPixelData(i)->GetCounts() >= fNumHitsMinimum ) {
	        fPointsId.push_back(i);
    	        fNumPoints++;
    	    }
    	}
    
        Msg(EsafMsg::Info) << "Num. of pixels with at least " << fNumHitsMinimum << " hits = " << fNumPoints << MsgDispatch;
    	
    } else if(! (gcm == NULL) ) {
        if ( gcm->GetObj("CluPixels") == NULL ) {
            Msg(EsafMsg::Warning) << "No data in event" << MsgDispatch;
	    return kFALSE;
        }
        fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
        if (fPointsId.size() == 0) {
            fNumPoints = 0;
	    Msg(EsafMsg::Warning) << "No pixels selected: terminated." << MsgDispatch;
            return kFALSE;
        } else {
            Msg(EsafMsg::Info) << "Pattern Recognition GtuClusteringModule data found with " << fPointsId.size()<< " points" << MsgDispatch;
            fNumPoints = fPointsId.size();
	}
    } else if ( !(moduledataHough->GetObj("Clusters") == NULL) ) {
        vector<EusoCluster*> cluHough;
        cluHough= *(vector<EusoCluster*>*) moduledataHough->GetObj("Clusters");
	fPointsId = cluHough[0]->GetPixelVector();
	if (fPointsId.size() == 0) {
            fNumPoints = 0;
            Msg(EsafMsg::Warning) << "No data in event" << MsgDispatch;
            return kFALSE;
        } else {
	    Msg(EsafMsg::Info) << "Pattern Recognition HoughModule data found with " << fPointsId.size() << " points" << MsgDispatch;
            fNumPoints = fPointsId.size();
	}
    }
	
    fNumHits = 0;
    for(Int_t i=0; i<fNumPoints; i++) {
        RecoPixelData *pix = fEv->GetRecoPixelData( fPointsId[i] );
        TVector3 pos(0,0,0);
        Double_t phitmp = pix->GetPhi()+PI;
        pos.SetMagThetaPhi( 1,pix->GetTheta(),phitmp );
        fData.fSpPos.push_back( pos );
        Int_t pixid = pix->GetPixelId();
        fData.fPixXYZ.push_back( runpars->PixelCenter(pixid) );
	fData.fSigmaPhi.push_back( pix->GetPhiSigma() );
	fData.fSigmaTheta.push_back( pix->GetThetaSigma() );
        fData.fTime.push_back( 0.01*fGtuLength*(Double_t)pix->GetGtu() ); // microseconds/100
        fData.fHits.push_back( pix->GetCounts() );
	fNumHits += pix->GetCounts();
    }
    fData.fNumPoints = fNumPoints;
    fData.fNumHits = fNumHits;
    fData.fGtuLength = fGtuLength;
	
    Msg(EsafMsg::Info) << "Processing " << fNumPoints << " pixel and " << fNumHits << " hits" << MsgDispatch;
	
    if(fData.fNumPoints < fNumPointsMin){
	fQuality = -1;
	Msg(EsafMsg::Warning) << "Not enough pixels selected ( minimum is "<< fNumPointsMin <<" ): terminate this event." << MsgDispatch;
        return kFALSE;
    } else {
	fQuality = 0;
    }
	
    fTrueTheta = fEv->GetHeader().GetTrueTheta();
    fTruePhi = fEv->GetHeader().GetTruePhi();
    fTrueDir.SetMagThetaPhi( 1., fTrueTheta, fTruePhi);
    Msg(EsafMsg::Info) << "TRUTH: Theta = " << fTrueTheta*DEG << " deg, Phi = " << fTruePhi*DEG << " deg, Energy = " << header.GetTrueEnergy() << " MeV." << MsgDispatch;
	
    TVector3 fTrueMaxOldsdr = fEv->GetHeader().GetTrueShowerMaxPos();
    Double_t xmax=fTrueMaxOldsdr.x()/km;
    Double_t ymax=fTrueMaxOldsdr.y()/km;
    Double_t zmax=fHISS-(fTrueMaxOldsdr.z()/km);
    fTrueMax.SetXYZ(xmax,ymax,zmax);
    Double_t phiMaxTmp=fTrueMax.Phi();
    fTrueMax.SetPhi(phiMaxTmp+PI);
    fTrueNorm = fTrueMax.Cross(fTrueDir); 
    if ( fTrueNorm.Z() > 0 ) fTrueNorm *= -1;
	
    Bool_t fShapeSelectionDone = kFALSE;
    if ( fDoHough ) {
	UseHoughandFindPlane();
	fShapeSelectionDone = kFALSE;
    } else if ( !fDoShapeSelection ) {
	FindPlane();
	fShapeSelectionDone = kFALSE;
    } else if ( fDoShapeSelection ) {
	UseShapeandFindPlane();
	if( fNumPointsSel < fNumPointsMin ){
	    Msg(EsafMsg::Info) << "Impossible to use selection by shape: only " << fNumPointsSel << " pixels selected!" <<MsgDispatch;
	    fShapeSelectionDone = kFALSE;
	    FindPlane();
	} else { fShapeSelectionDone = kFALSE; }
    } 
	
    Double_t fangleDirNorm = fTrueDir.Angle(fNorm);
    fDeltaTDP = fTrueNorm.Angle(fNorm);
	
    Msg(EsafMsg::Info) << "Found TrackDetectorPlane with error = " << fDeltaTDP*DEG << " deg, fangleDirNorm = " << fangleDirNorm*DEG << " deg" << MsgDispatch;
	
    fW = (TVector3(0,0,1).Cross(fNorm)).Unit();		
    fU = (fNorm.Cross(fW)).Unit();					
	
	
    if( fOptionSelectionHough || (fUseShapeSelectioninModules &&  fShapeSelectionDone ) ){
	fNumPoints =  fNumPointsSel;
	fNumHits =  fNumHitsSel;
	fData.fNumPoints = fNumPoints;
	fData.fNumHits = fNumHits;
	fData.fSpPos.clear();
	fData.fTime.clear();
	fData.fHits.clear();
	fData.fSigmaTheta.clear();
	fData.fSigmaPhi.clear();
	for(Int_t i=0;i<fNumPointsSel;i++){
	    fData.fSpPos.push_back(fData.fSpPosSel[i]);
	    fData.fTime.push_back(fData.fTimeSel[i]);
	    fData.fHits.push_back(fData.fHitsSel[i]);
	    fData.fSigmaTheta.push_back(fData.fSigmaThetaSel[i]);
	    fData.fSigmaPhi.push_back(fData.fSigmaPhiSel[i]);
	}
    }

    switch (fMethodIdentifier) {
        case 1: AA1(); break;
        case 2: NE1(); break;
        case 3: AA2(); break;
        case 4: NE2(); break;
        case 5: AE1(); break;
        case 6: all(); break;
        default: Msg(EsafMsg::Panic) << "Wrong fMethod config value in TrackDirection2Module.cfg" << MsgDispatch; break;
    } 
	
    fData.Clear();
    
    //save data
    MyData()->Add("NormX",fNorm.X());
    MyData()->Add("NormY",fNorm.Y());
    MyData()->Add("NormZ",fNorm.Z());
    MyData()->Add("WAxisX",fW.X());
    MyData()->Add("WAxisY",fW.Y());
    MyData()->Add("WAxisZ",fW.Z());
    MyData()->Add("Theta",fTHETAreco);
    MyData()->Add("Phi",fPHIreco);
	
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t TrackDirection2Module::PostProcess() {
    //
    // Post-processing method
    //
    
    if( fQuality == 0) {
	fRecoEventsCounter++;
	vecDeltaTDP.push_back(fDeltaTDP*DEG);
	if( fMethodIdentifier == 6 ) {	
	    vecDeltaEASDirAA1.push_back(fDeltaEASDirAA1*DEG);
	    vecDeltaEASDirAA2.push_back(fDeltaEASDirAA2*DEG);
	    vecDeltaEASDirAE1.push_back(fDeltaEASDirAE1*DEG);
	    vecDeltaEASDirNE1.push_back(fDeltaEASDirNE1*DEG);
	    vecDeltaEASDirNE2.push_back(fDeltaEASDirNE2*DEG);
	} else {	
	    vecDeltaEASDir.push_back(fDeltaEASDir*DEG);
	}
    }		
			
    return kTRUE;

}

//_____________________________________________________________________________
 Bool_t TrackDirection2Module::Done() { 
    //
    // Module done method. Do some statistics about the renconstructed events
    //

    Msg(EsafMsg::Info) << "Number of reconstructed events = " << fRecoEventsCounter << MsgDispatch;
    if (fRecoEventsCounter < 2) return kTRUE;
	
    gDirectory->cd("/");
    string pathname = "TD2Resolution";
    TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
    path->cd();
	
    string aa1 = "AA1";
    string aa2 = "AA2";	
    string ae1 = "AE1";
    string ne1 = "NE1";
    string ne2 = "NE2";
    string tdp = "TDP";
	
    DoStat(vecDeltaTDP,tdp);
    if ( fMethodIdentifier == 6 ){
	DoStat(vecDeltaEASDirAA1,aa1);
	DoStat(vecDeltaEASDirAA2,aa2);
	DoStat(vecDeltaEASDirAE1,ae1);
	DoStat(vecDeltaEASDirNE1,ne1);
	DoStat(vecDeltaEASDirNE2,ne2);
    } else {
	if ( fMethodIdentifier == 1 ) DoStat(vecDeltaEASDir,aa1);
	if ( fMethodIdentifier == 2 ) DoStat(vecDeltaEASDir,ne1);
	if ( fMethodIdentifier == 3 ) DoStat(vecDeltaEASDir,aa2);
	if ( fMethodIdentifier == 4 ) DoStat(vecDeltaEASDir,ne2);
	if ( fMethodIdentifier == 5 ) DoStat(vecDeltaEASDir,ae1);
    }
	
    vecDeltaEASDir.clear();
    vecDeltaEASDirAA1.clear();
    vecDeltaEASDirAA2.clear();
    vecDeltaEASDirAE1.clear();
    vecDeltaEASDirNE1.clear();
    vecDeltaEASDirNE2.clear();
    vecDeltaTDP.clear();
    fRecoEventsCounter = 0;

    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
    return kTRUE;

}

//_____________________________________________________________________________
 void TrackDirection2Module::DoStat(vector<Double_t> err, string namemethod){
    //
    // Statistics of reconstructed events.
    //

    Int_t best(0);
    Int_t med(0);
    Int_t worst(0);
    Int_t numbins=400;
    Int_t maxDeltaEAS=20;
    TH1F *hR=new TH1F("hR","hR",numbins,0,maxDeltaEAS);
	
    for( Int_t i=0; i<fRecoEventsCounter; i++ ) {
	hR->Fill(err[i]);
	if ( err[i] < fStat1 ) best++;
	if ( err[i] >= fStat1 && err[i] <= fStat2 ) med++;
	if ( err[i] > fStat2 ) worst++;
    }
	
    Bool_t foundresolution = kFALSE;
    Stat_t under68(0);
    Float_t frac(0);
    Double_t resolution(0);
    for(Int_t i=0;i<numbins;i++){
	under68 += hR->GetBinContent(i);
	frac=100*(Float_t)under68/((Float_t)fRecoEventsCounter);
	if(frac>68){
	    resolution=i*maxDeltaEAS/((Float_t)numbins);
	    Msg(EsafMsg::Info)<<"Eas-Dir reconstruction ("<<namemethod.c_str()<<"):tResolution(68%)="<<resolution<<" deg"<<MsgDispatch;
	    foundresolution = kTRUE;
	    break;
	}
    }
    if (!foundresolution) 
        Msg(EsafMsg::Info)<<"Eas-Dir reconstruction (" << namemethod.c_str() <<"):tResolution(68%) not found, too bad!!" << MsgDispatch;
	
    string titleh = namemethod.c_str();
    titleh += Form("  Resolution(68%)=%g deg", resolution);
    hR->SetTitle(titleh.c_str());
    hR->GetXaxis()->SetTitle("#Psi (deg)");
    hR->Write();
    delete hR;
	
    Msg(EsafMsg::Info) << "tEvents with error under " << fStat1<< " deg = " <<best<< " ( " << ((Float_t)best/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
    Msg(EsafMsg::Info) << "tEvents with error between "<<fStat1<<" and "<<fStat2<<" deg = " <<med<< " ( " << ((Float_t)med/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
    Msg(EsafMsg::Info) << "tEvents with error above "<<fStat2<<" deg = " <<worst<< " ( " << ((Float_t)worst/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
	
    return;
}

//_____________________________________________________________________________
 void TrackDirection2Module::UserMemoryClean()  {
    //
    // User memory clean
    //
}

//_____________________________________________________________________________
 Bool_t TrackDirection2Module::SaveRootData(RecoRootEvent *fRecoRootEvent) {
    //
    // Save data in the reco rootfile
    //
    
    fRecoRootEvent->GetRecoTrackDirection2().SetErrorTDP(fDeltaTDP*DEG);
    fRecoRootEvent->GetRecoTrackDirection2().SetQuality(fQuality); // -1 if reconstruction is not possible in this context, otherwise 0 
    
    if(fMethodIdentifier == 6){
    	fRecoRootEvent->GetRecoTrackDirection2().SetAA1(fTHETArecoAA1*DEG, fPHIrecoAA1*DEG, fDeltaEASDirAA1*DEG);
    	fRecoRootEvent->GetRecoTrackDirection2().SetAA2(fTHETArecoAA2*DEG, fPHIrecoAA2*DEG, fDeltaEASDirAA2*DEG);
    	fRecoRootEvent->GetRecoTrackDirection2().SetNE1(fTHETArecoNE1*DEG, fPHIrecoNE1*DEG, fDeltaEASDirNE1*DEG);
    	fRecoRootEvent->GetRecoTrackDirection2().SetNE2(fTHETArecoNE2*DEG, fPHIrecoNE2*DEG, fDeltaEASDirNE2*DEG);
    	fRecoRootEvent->GetRecoTrackDirection2().SetAE1(fTHETArecoAE1*DEG, fPHIrecoAE1*DEG, fDeltaEASDirAE1*DEG);
    }else{
    	fRecoRootEvent->GetRecoTrackDirection2().SetTheta(fTHETAreco*DEG);
    	fRecoRootEvent->GetRecoTrackDirection2().SetErrorTheta(fDeltaTheta*DEG);
    	fRecoRootEvent->GetRecoTrackDirection2().SetPhi(fPHIreco*DEG); 
    	fRecoRootEvent->GetRecoTrackDirection2().SetErrorPhi(fDeltaPhi*DEG);
    	fRecoRootEvent->GetRecoTrackDirection2().SetErrorDirection(fDeltaEASDir*DEG);
    }
    
    return kTRUE;
}

//______________________________________________________________________________
 void TrackDirection2Module::UseHoughandFindPlane(){
    //
    // Find track-detector plane (TDP) using Hough Transform
    //
	
    Msg(EsafMsg::Info) << "UseHoughandFindPlane()" << MsgDispatch;	
	
    vector<Double_t> x;
    vector<Double_t> y;
    vector<Double_t> t;
    vector<Int_t> c;
    for( Int_t i=0; i<fNumPoints; i++ ) {
	TVector3 dummy = fData.fSpPos[i];
	t.push_back(fData.fTime[i]);
  	x.push_back(dummy.x()/dummy.z());
        y.push_back(dummy.y()/dummy.z());
        c.push_back(fData.fHits[i]);
    }
    
    // errors for the Hough fit
    Float_t ex = 0.001;
    Float_t ey = 0.001;
    Float_t et = 2.5*0.01;
    
    // use Hough fit to select track points
    // x-t Hough fit
    HoughFit *xtFit = new HoughFit( t, x, c, fOptionSelectionHough,et,ex,fPointsId);
    Double_t vX = xtFit->GetSlope();
    Double_t wX = xtFit->GetOffset();
    vector<Int_t> idSelX = xtFit->GetIdSel();
    delete xtFit;
    
    // y-t Hough fit
    HoughFit *ytFit = new HoughFit( t, y, c, fOptionSelectionHough,et,ey,fPointsId);
    Double_t vY = ytFit->GetSlope();
    Double_t wY = ytFit->GetOffset();
    vector<Int_t> idSelY = ytFit->GetIdSel();
    delete ytFit;
	
    // compute the normal vector to the TDP
    fNorm.SetXYZ(vY,-vX,wY*vX-wX*vY); 
    fNorm.SetMag(1.);
    if ( fNorm.Z() > 0 ) fNorm *= -1;
		
    //use the following only to debug (to see how the projection on plane Z=1 versus time appear) 
    if( fDoGraphUseShape ){
        Double_t xprog[fNumPoints], yprog[fNumPoints], tprog[fNumPoints];
        Double_t xline[fNumPoints],yline[fNumPoints];
	for( Int_t i=0; i<fNumPoints; i++ ) {
	    TVector3 dummy = fData.fSpPos[i];
	    tprog[i] = fData.fTime[i];
	    xline[i] = vX*fData.fTime[i] + wX;
	    yline[i]=vY*fData.fTime[i]+wY;
	    xprog[i] = dummy.x()/dummy.z();
	    yprog[i] = dummy.y()/dummy.z();
	}
	TGraph *gxypro = new TGraph(fNumPoints,yprog,xprog);
	gxypro->SetNameTitle("gxypro","xprog vs yprog");gxypro->SetMarkerSize(.5);gxypro->SetMarkerStyle(23);
	TGraph *gxpro = new TGraph(fNumPoints,tprog,xprog);
	gxpro->SetNameTitle("gxpro","xprog vs time");gxpro->SetMarkerSize(.5);gxpro->SetMarkerStyle(23);
	TGraph *gypro = new TGraph(fNumPoints,tprog,yprog);
	gypro->SetNameTitle("gypro","yprog vs time");gypro->SetMarkerSize(.5);gypro->SetMarkerStyle(23);
	TGraph *gxproretta = new TGraph(fNumPoints,tprog,xline);
	gxproretta->SetMarkerColor(4);gxproretta->SetLineColor(4);
	TGraph *gyproretta = new TGraph(fNumPoints,tprog,yline);
	gyproretta->SetMarkerColor(4);gyproretta->SetLineColor(4);
	TCanvas *cUSFP=new TCanvas("cUSFP","cUSFP",1);
	cUSFP->Divide(3,1);
	cUSFP->cd(1);gxpro->Draw("AP");gxproretta->Draw("PLsame");
	cUSFP->cd(2);gypro->Draw("AP");gyproretta->Draw("PLsame");
	cUSFP->cd(3);gxypro->Draw("AP");
	cUSFP->Write();
    }
	
    //cout<<"fPointsId.size()="<<fPointsId.size()<<"tidSelX.size()="<<idSelX.size()<<"tidSelY.size()="<<idSelY.size()<<endl;
    vector<Int_t> idSelectedMerged(idSelY.size()+idSelX.size());
    merge(idSelX.begin(), idSelX.end(),idSelY.begin(), idSelY.end(),idSelectedMerged.begin());
    vector<Int_t>::iterator different = unique(idSelectedMerged.begin(),idSelectedMerged.end());
    idSelectedMerged.erase(different, idSelectedMerged.end());
	
    fNumPointsSel = idSelectedMerged.size();
    Msg(EsafMsg::Info) << "Selected " << fNumPointsSel << " points over " << fPointsId.size() << MsgDispatch;
    
    // if there are > 10 points selected recompute the TDP more precisely
    fNumHitsSel = 0;
    if( fNumPointsSel > 10 ){
	vector<Int_t> csel;
	vector<Double_t> xsel;
	vector<Double_t> ysel;
	vector<Double_t> errxsel;
	vector<Double_t> errysel;
	vector<Double_t> tsel;
	vector<Double_t> xsel1;
	vector<Double_t> ysel1;
        vector<Double_t> tsel1;
		
	for(Int_t i=0; i<fNumPointsSel; i++) {
	    RecoPixelData *pix = fEv->GetRecoPixelData( idSelectedMerged[i] );
	    TVector3 pos(0,0,0);
    	    Double_t phitmp = pix->GetPhi() + PI;
    	    pos.SetMagThetaPhi( 1, pix->GetTheta(), phitmp );
    	    xsel.push_back( pos.x()/pos.z() );
    	    ysel.push_back( pos.y()/pos.z() );
    	    tsel.push_back( 0.01*(Double_t)2.500*pix->GetGtu() );
    	    csel.push_back( pix->GetCounts() );
    	    fNumHitsSel += csel[i];
    	    Double_t errphi = pix->GetPhiSigma();
    	    Double_t errtheta = pix->GetThetaSigma();
    	    fData.fSpPosSel.push_back(pos); 
	    fData.fTimeSel.push_back(tsel[i]); 
	    fData.fHitsSel.push_back(csel[i]);
	    fData.fSigmaThetaSel.push_back(errtheta); 
	    fData.fSigmaPhiSel.push_back(errphi); 
            Double_t dx = DeltaX(pos.Theta(), pos.Phi(), errtheta, errphi);
            Double_t dy = DeltaY(pos.Theta(), pos.Phi(), errtheta, errphi);
	    errxsel.push_back(dx/csel[i]);
            errysel.push_back(dy/csel[i]);
            // add hits for median fit
            for (Int_t j(0); j<pix->GetCounts(); j++) {
    	        xsel1.push_back( pos.x()/pos.z() );
    	        ysel1.push_back( pos.y()/pos.z() );
    	        tsel1.push_back( 0.01*(Double_t)2.500*pix->GetGtu() );
            }
	}
	
        TVector3 fNormA, fNormB, fNormC;
        Double_t vX2(0), vY2(0), wX2(0), wY2(0);
        if ( fDebugInfo || useHFplane ) { // recalculate normal to TDP using Hough fit
	    HoughFit *xhf1 = new HoughFit( tsel, xsel, csel, 0,et,ex);
    	    vX2 = xhf1->GetSlope();
            wX2 = xhf1->GetOffset();
    	    delete xhf1;
	    HoughFit *yhf1 = new HoughFit( tsel, ysel, csel, 0,et,ey);
    	    vY2 = yhf1->GetSlope();
    	    wY2 = yhf1->GetOffset();
	    delete yhf1;
        
	    fNormA.SetXYZ(vY2,-vX2,wY2*vX2-wX2*vY2); 
   	    fNormA.SetMag(1.);
            if ( fNormA.Z() > 0 ) fNormA *= -1;
            Double_t fDeltaTDPhh = fTrueNorm.Angle(fNormA);
	    Msg(EsafMsg::Info) << "TDP error (hough selection + hough fit) = " << fDeltaTDPhh*DEG << MsgDispatch;
	}
        
        if ( fDebugInfo || useLFplane ) { // recalculate normal to TDP using least squares fit
	    LeastSquaresFit *xlf1 = new LeastSquaresFit( fNumPointsSel, tsel, xsel, errxsel );
    	    Double_t vX2ls = xlf1->GetSlope();
    	    Double_t wX2ls = xlf1->GetOffset();
	    delete xlf1;
	    LeastSquaresFit *ylf1 = new LeastSquaresFit( fNumPointsSel, tsel, ysel, errysel );
    	    Double_t vY2ls = ylf1->GetSlope();
    	    Double_t wY2ls = ylf1->GetOffset();
	    delete ylf1;
	
	    fNormB.SetXYZ(vY2ls,-vX2ls,wY2ls*vX2ls-wX2ls*vY2ls); 
   	    fNormB.SetMag(1.);	if ( fNormB.Z() > 0 ) fNormB *= -1;
	    Double_t fDeltaTDPhls = fTrueNorm.Angle(fNormB);
	    Msg(EsafMsg::Info) << "TDP error (hough selection + linear fit) = " << fDeltaTDPhls*DEG << MsgDispatch;
        }
	    
        if ( fDebugInfo || useHFplane ) { // recalculate normal to TDP using median fit
	    MedianFit *xmf1 = new MedianFit( fNumHitsSel, tsel1, xsel1 ); // CHANGED
    	    Double_t vX2m = xmf1->GetSlope();
    	    Double_t wX2m = xmf1->GetOffset();
	    delete xmf1;
	    MedianFit *ymf1 = new MedianFit( fNumHitsSel, tsel1, ysel1 ); // CHANGED
    	    Double_t vY2m = ymf1->GetSlope();
    	    Double_t wY2m = ymf1->GetOffset();
	    delete ymf1;
        
	    fNormC.SetXYZ(vY2m,-vX2m,wY2m*vX2m-wX2m*vY2m); 
 	    fNormC.SetMag(1.); 		
            if ( fNormC.Z() > 0 ) fNormC *= -1;
	    Double_t fDeltaTDPhm = fTrueNorm.Angle(fNormC);
	    Msg(EsafMsg::Info) << "TDP error (hough selection + median fit) = " << fDeltaTDPhm*DEG << MsgDispatch;
        }
	
        if ( useLFplane ) fNorm = fNormB;
    	else if ( useMFplane ) fNorm = fNormC;
        else if ( useHFplane ) fNorm = fNormA;
        fNorm.SetMag(1.);	
        if ( fNorm.Z() > 0 ) fNorm *= -1;
	
	//use the following only to debug (to see how the projection on plane Z=1 versus time appear) 
	if( fDoGraphUseShape ){
	    Double_t xprog2[fNumPointsSel],yprog2[fNumPointsSel],tprog2[fNumPointsSel],xline2[fNumPointsSel],yline2[fNumPointsSel];
	    for( Int_t i=0; i<fNumPointsSel; i++ ) {
		tprog2[i] = tsel[i];
		xline2[i] = vX2*tsel[i]+wX2;
		yline2[i] = vY2*tsel[i]+wY2;
		xprog2[i] = xsel[i];
		yprog2[i] = ysel[i];
	    }
	    TGraph *gxypro2 = new TGraph(fNumPointsSel,yprog2,xprog2);
	    gxypro2->SetNameTitle("gxypro2","xprog2 vs yprog2");gxypro2->SetMarkerSize(.5);gxypro2->SetMarkerStyle(23);
	    TGraph *gxpro2 = new TGraph(fNumPointsSel,tprog2,xprog2);
	    gxpro2->SetNameTitle("gxpro2","xprog2 vs time2");gxpro2->SetMarkerSize(.5);gxpro2->SetMarkerStyle(23);
	    TGraph *gypro2 = new TGraph(fNumPointsSel,tprog2,yprog2);
	    gypro2->SetNameTitle("gypro2","yprog2 vs time2");gypro2->SetMarkerSize(.5);gypro2->SetMarkerStyle(23);
	    TGraph *gxproretta2 = new TGraph(fNumPointsSel,tprog2,xline2);
	    gxproretta2->SetMarkerColor(4);gxproretta2->SetLineColor(4);
	    TGraph *gyproretta2 = new TGraph(fNumPointsSel,tprog2,yline2);
	    gyproretta2->SetMarkerColor(4);gyproretta2->SetLineColor(4);
	    TCanvas *cUSFP2=new TCanvas("cUSFP2","cUSFP2",1);
	    cUSFP2->Divide(3,1);
	    cUSFP2->cd(1);gxpro2->Draw("AP");gxproretta2->Draw("PLsame");
	    cUSFP2->cd(2);gypro2->Draw("AP");gyproretta2->Draw("PLsame");
	    cUSFP2->cd(3);gxypro2->Draw("AP");
	    cUSFP2->Write();
	}
    }
}

//______________________________________________________________________________
 void TrackDirection2Module::FindPlane(){
    //
    // Find the track-detector plane (TDP)
    //
    
    Msg(EsafMsg::Info) << "FindPlane()" << MsgDispatch;	

    vector<Double_t> xpro;
    vector<Double_t> ypro;
    vector<Double_t> errxpro;
    vector<Double_t> errypro;
    vector<Double_t> tpro; 
    vector<Int_t> cpro;
    vector<Double_t> xpro1;
    vector<Double_t> ypro1;
    vector<Double_t> tpro1; 
	
    for( Int_t i=0; i<fNumPoints; i++ ) {
	TVector3 dummy = fData.fSpPos[i];
        Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
        Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
	xpro.push_back(dummy.X()/dummy.Z());
        ypro.push_back(dummy.Y()/dummy.Z());
        cpro.push_back(fData.fHits[i]);
        tpro.push_back(fData.fTime[i]);
	errxpro.push_back(dx/fData.fHits[i]);
        errypro.push_back(dy/fData.fHits[i]);
        // add hits for median fit
        for( Int_t j(0); j<fData.fHits[i]; j++ ) {
 	    xpro1.push_back(dummy.X()/dummy.Z());
            ypro1.push_back(dummy.Y()/dummy.Z());
            tpro1.push_back(fData.fTime[i]);           
        }
    }

    TVector3 fNormA, fNormB, fNormC;
    if ( fDebugInfo || useHFplane ) { // calculate normal to TDP using Hough fit
        // errors for the Hough fit
        Float_t ex = 0.001;
        Float_t ey = 0.001;
        Float_t et = 2.5*0.01;
        HoughFit *xhf1 = new HoughFit( tpro, xpro, cpro, 0,et,ex);
    	Double_t vX2 = xhf1->GetSlope();
        Double_t wX2 = xhf1->GetOffset();
    	delete xhf1;
	HoughFit *yhf1 = new HoughFit( tpro, ypro, cpro, 0,et,ey);
    	Double_t vY2 = yhf1->GetSlope();
    	Double_t wY2 = yhf1->GetOffset();
	delete yhf1;
        
	fNormA.SetXYZ(vY2,-vX2,wY2*vX2-wX2*vY2); 
   	fNormA.SetMag(1.);
        if ( fNormA.Z() > 0 ) fNormA *= -1;
        Double_t fDeltaTDPhh = fTrueNorm.Angle(fNormA);
	Msg(EsafMsg::Info) << "TDP error (hough fit) = " << fDeltaTDPhh*DEG << MsgDispatch;
    }
        
    if ( fDebugInfo || useLFplane ) { // calculate normal to TDP using least squares fit
	LeastSquaresFit *xlf1 = new LeastSquaresFit( fNumPoints, tpro, xpro, errxpro );
    	Double_t vX2ls = xlf1->GetSlope();
    	Double_t wX2ls = xlf1->GetOffset();
	delete xlf1;
	LeastSquaresFit *ylf1 = new LeastSquaresFit( fNumPoints, tpro, ypro, errypro );
    	Double_t vY2ls = ylf1->GetSlope();
    	Double_t wY2ls = ylf1->GetOffset();
	delete ylf1;
	
	fNormB.SetXYZ(vY2ls,-vX2ls,wY2ls*vX2ls-wX2ls*vY2ls); 
   	fNormB.SetMag(1.);	if ( fNormB.Z() > 0 ) fNormB *= -1;
	Double_t fDeltaTDPhls = fTrueNorm.Angle(fNormB);
	Msg(EsafMsg::Info) << "TDP error (linear fit) = " << fDeltaTDPhls*DEG << MsgDispatch;
    }
	    
    if ( fDebugInfo || useHFplane ) { // recalculate normal to TDP using median fit
	MedianFit *xmf1 = new MedianFit( fNumHits, tpro1, xpro1 );
    	Double_t vX2m = xmf1->GetSlope();
    	Double_t wX2m = xmf1->GetOffset();
	delete xmf1;
	MedianFit *ymf1 = new MedianFit( fNumHits, tpro1, ypro1 );
    	Double_t vY2m = ymf1->GetSlope();
    	Double_t wY2m = ymf1->GetOffset();
	delete ymf1;
        
	fNormC.SetXYZ(vY2m,-vX2m,wY2m*vX2m-wX2m*vY2m); 
 	fNormC.SetMag(1.); 		
        if ( fNormC.Z() > 0 ) fNormC *= -1;
	Double_t fDeltaTDPhm = fTrueNorm.Angle(fNormC);
	Msg(EsafMsg::Info) << "TDP error (median fit) = " << fDeltaTDPhm*DEG << MsgDispatch;
    }
	
    if ( useLFplane ) fNorm = fNormB;
    else if ( useMFplane ) fNorm = fNormC;
    else if ( useHFplane ) fNorm = fNormA;
    fNorm.SetMag(1.);
    if ( fNorm.Z() > 0 ) fNorm *= -1;

    cpro.clear();
    xpro.clear();
    ypro.clear();
    tpro.clear();
    errxpro.clear();
    errypro.clear();
}

//______________________________________________________________________________
 void TrackDirection2Module::UseShapeandFindPlane(){
    //
    // Find the track-detector plane (TDP) using the shape method
    //
    
    vector<Int_t> cpro;
    vector<Double_t> xpro;
    vector<Double_t> ypro;
    vector<Double_t> errxpro;
    vector<Double_t> errypro;
    vector<Double_t> tpro; 
    vector<Double_t> xpro1;
    vector<Double_t> ypro1;
    vector<Double_t> tpro1;
    
    for( Int_t i=0; i<fNumPoints; i++ ) {
	TVector3 dummy = fData.fSpPos[i];
        Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
        Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
	xpro.push_back(dummy.X()/dummy.Z());
        ypro.push_back(dummy.Y()/dummy.Z());
        tpro.push_back(fData.fTime[i]);
        cpro.push_back(fData.fHits[i]);
	errxpro.push_back(dx/fData.fHits[i]);
        errypro.push_back(dy/fData.fHits[i]);
        // add hits for median fit
        for(Int_t j(0); j<fData.fHits[i]; i++) {
	    xpro1.push_back(dummy.X()/dummy.Z());
            ypro1.push_back(dummy.Y()/dummy.Z());
            tpro1.push_back(fData.fTime[i]);
        }
    }
    
    // FIRST STEP: median fit of points
    MedianFit *xMf1 = new MedianFit(fNumHits, tpro1, xpro1); //CHANGED
    Double_t sX = xMf1->GetSlope();
    Double_t oX = xMf1->GetOffset();
    Double_t aX = xMf1->GetAbsoluteDeviation();
    delete xMf1;
    MedianFit *yMf1 = new MedianFit(fNumHits, tpro1, ypro1); //CHANGED
    Double_t sY = yMf1->GetSlope();
    Double_t oY = yMf1->GetOffset();
    Double_t aY = yMf1->GetAbsoluteDeviation();
    delete yMf1;
	
    //use the following only to debug (to see how the projection on plane Z=1 versus time appear) 
    if( fDoGraphUseShape ){
	Double_t xprog[fNumPoints], yprog[fNumPoints], tprog[fNumPoints];
        Double_t xline[fNumPoints], yline[fNumPoints];
	Double_t xlineup[fNumPoints], ylineup[fNumPoints];
	Double_t xlinedown[fNumPoints], ylinedown[fNumPoints];
	for( Int_t i=0; i<fNumPoints; i++ ) {
	    TVector3 dummy = fData.fSpPos[i];
	    tprog[i] = fData.fTime[i];
	    xline[i] = sX * fData.fTime[i] + oX;
	    xlineup[i] = xline[i] + fMulti1 * aX;
	    xlinedown[i] = xline[i] - fMulti1 * aX;
	    yline[i] = sY * fData.fTime[i] + oY;
	    ylineup[i] = yline[i] + fMulti1 * aY;
	    ylinedown[i] = yline[i] - fMulti1 * aY;
	    xprog[i] = dummy.X()/dummy.Z();
	    yprog[i] = dummy.Y()/dummy.Z();
	}
	TGraph *gxypro = new TGraph(fNumPoints,yprog,xprog);
	gxypro->SetNameTitle("gxypro","xprog vs yprog");gxypro->SetMarkerSize(.5);gxypro->SetMarkerStyle(23);
	TGraph *gxpro = new TGraph(fNumPoints,tprog,xprog);
	gxpro->SetNameTitle("gxpro","xprog vs time");gxpro->SetMarkerSize(.5);gxpro->SetMarkerStyle(23);
	TGraph *gypro = new TGraph(fNumPoints,tprog,yprog);
	gypro->SetNameTitle("gypro","yprog vs time");
	gypro->SetMarkerSize(.5);
	gypro->SetMarkerStyle(23);
	TGraph *gxproretta = new TGraph(fNumPoints,tprog,xline);
	gxproretta->SetMarkerColor(4);gxproretta->SetLineColor(4);
	TGraph *gyproretta = new TGraph(fNumPoints,tprog,yline);
	gyproretta->SetMarkerColor(4);gyproretta->SetLineColor(4);
	TGraph *gxprorettaup = new TGraph(fNumPoints,tprog,xlineup);
	gxprorettaup->SetMarkerColor(2);gxprorettaup->SetLineColor(2);
	TGraph *gyprorettaup = new TGraph(fNumPoints,tprog,ylineup);
	gyprorettaup->SetMarkerColor(2);gyprorettaup->SetLineColor(2);
	TGraph *gxprorettadown = new TGraph(fNumPoints,tprog,xlinedown);
	gxprorettadown->SetMarkerColor(2);gxprorettadown->SetLineColor(2);
	TGraph *gyprorettadown = new TGraph(fNumPoints,tprog,ylinedown);
	gyprorettadown->SetMarkerColor(2);gyprorettadown->SetLineColor(2);
	TCanvas *cUSFP=new TCanvas("cUSFP","cUSFP",1);
	cUSFP->Divide(3,1);
	cUSFP->cd(1);gxpro->Draw("AP");gxproretta->Draw("PLsame");gxprorettaup->Draw("PLsame");gxprorettadown->Draw("PLsame");
	cUSFP->cd(2);gypro->Draw("AP");gyproretta->Draw("PLsame");gyprorettaup->Draw("PLsame");gyprorettadown->Draw("PLsame");
	cUSFP->cd(3);gxypro->Draw("AP");
	cUSFP->Write();
    }

    // make a shape selection of track points
    vector<Double_t> xprosel;
    vector<Double_t> yprosel;
    vector<Double_t> errxprosel;
    vector<Double_t> erryprosel;
    vector<Double_t> tprosel;
    vector<Double_t> xprosel1;
    vector<Double_t> yprosel1;
    vector<Double_t> tprosel1;
	
    vector<TVector3> fSpPosSeltmp;
    vector<Double_t> fTimeSeltmp;
    vector<Int_t> fHitsSeltmp;
    vector<Double_t> fSigmaThetaSeltmp;
    vector<Double_t> fSigmaPhiSeltmp;
	
    Int_t fNumPointsSeltmp = 0;
    Int_t fNumHitsSeltmp = 0;
    for( Int_t i=0; i<fNumPoints; i++ ) {
	TVector3 dummy = fData.fSpPos[i];
	Double_t xprotmp=dummy.X()/dummy.Z();
	Double_t yprotmp=dummy.Y()/dummy.Z();
	Double_t tprotmp=fData.fTime[i];
        // shape selection with multiplier fMulti1
	if( Abs(xprotmp-sX*tprotmp-oX) < fMulti1*aX || Abs(yprotmp-sY*tprotmp-oY) < fMulti1*aY) {
            Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
            Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
	    fSpPosSeltmp.push_back(dummy); 
	    fTimeSeltmp.push_back(fData.fTime[i]); 
	    fHitsSeltmp.push_back(fData.fHits[i]);
	    fSigmaThetaSeltmp.push_back(fData.fSigmaTheta[i]); 
	    fSigmaPhiSeltmp.push_back(fData.fSigmaPhi[i]); 
	    xprosel.push_back(xprotmp);
	    yprosel.push_back(yprotmp);
	    tprosel.push_back(tprotmp);
	    errxprosel.push_back(dx/fData.fHits[i]);
            erryprosel.push_back(dy/fData.fHits[i]);
	    fNumPointsSeltmp++;
	    fNumHitsSeltmp += fData.fHits[i];
            //add hits for median fit
            for( Int_t j(0); j<fData.fHits[i]; j++) {
    	        xprosel1.push_back(xprotmp);
    	        yprosel1.push_back(yprotmp);
	        tprosel1.push_back(tprotmp);
            }
	} 
    }
    Msg(EsafMsg::Info) << "UseShapeandFindPlane(), first step: selected " << fNumPointsSeltmp << " from " << fNumPoints << " pixels" << MsgDispatch;
	
    if ( fNumPointsSeltmp < fNumPointsMin ){ // return if not enough points
	fNumPointsSel = fNumPointsSeltmp;
	return;
    }

    // SECOND STEP: new median fit of selected points
    MedianFit *xMf2 = new MedianFit(fNumHitsSeltmp, tprosel1, xprosel1); //CHANGED
    Double_t sX2 = xMf2->GetSlope();
    Double_t oX2 = xMf2->GetOffset();
    Double_t aX2 = xMf2->GetAbsoluteDeviation();
    delete xMf2;
    MedianFit *yMf2 = new MedianFit(fNumHitsSeltmp, tprosel1, yprosel1); //CHANGED
    Double_t sY2 = yMf2->GetSlope();
    Double_t oY2 = yMf2->GetOffset();
    Double_t aY2 = yMf2->GetAbsoluteDeviation();
    delete yMf2;

    // new shape selection of points from median fit results
    vector<Double_t> xprosel2;
    vector<Double_t> yprosel2;
    vector<Double_t> errxprosel2;
    vector<Double_t> erryprosel2;
    vector<Double_t> tprosel2;
    vector<Int_t> cprosel2;
    vector<Double_t> xprosel3;
    vector<Double_t> yprosel3;
    vector<Double_t> tprosel3;
	
    fNumPointsSel = 0;
    fNumHitsSel = 0;
	
    for( Int_t i=0; i<fNumPointsSeltmp; i++ ) {
	TVector3 dummy = fSpPosSeltmp[i];
	Double_t xprotmp = dummy.X()/dummy.Z();
	Double_t yprotmp = dummy.Y()/dummy.Z();
	Double_t tprotmp = fTimeSeltmp[i];
        // shape selction with multiplier fMulti2
	if( Abs(xprotmp-sX2*tprotmp-oX2) < fMulti2*aX2 || TMath::Abs(yprotmp-sY2*tprotmp-oY2) < fMulti2*aY2 ){
            Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
            Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
	    fData.fSpPosSel.push_back(dummy); 
	    fData.fTimeSel.push_back(fTimeSeltmp[i]); 
	    fData.fHitsSel.push_back(fHitsSeltmp[i]);
	    fData.fSigmaThetaSel.push_back(fSigmaThetaSeltmp[i]); 
	    fData.fSigmaPhiSel.push_back(fSigmaPhiSeltmp[i]); 
	    cprosel2.push_back(fHitsSeltmp[i]);
	    xprosel2.push_back(xprotmp);
	    yprosel2.push_back(yprotmp);
	    tprosel2.push_back(tprotmp);
	    errxprosel2.push_back(dx/fHitsSeltmp[i]);
            erryprosel2.push_back(dy/fHitsSeltmp[i]);
	    fNumPointsSel++;
	    fNumHitsSel += fHitsSeltmp[i];
            //add hits for median fit
            for (Int_t j(0); j<fHitsSeltmp[i]; j++) {
    	        xprosel3.push_back(xprotmp);
	        yprosel3.push_back(yprotmp);
	        tprosel3.push_back(tprotmp);
            }
	} 
    }
    Msg(EsafMsg::Info) << "UseShapeandFindPlane(), second step: selected " << fNumPointsSel << " from " << fNumPointsSeltmp << " pixels" << MsgDispatch;

    if ( fNumPointsSel < fNumPointsMin ) return; // return if not enough points

    Double_t vX(0), wX(0), vY(0), wY(0);
    // fit of selected points
    if ( useLFplane ) {
        LeastSquaresFit *xtFit = new LeastSquaresFit(fNumPointsSel, tprosel2, xprosel2, errxprosel2);
        vX = xtFit->GetSlope();
        wX = xtFit->GetOffset();
        delete xtFit;
        LeastSquaresFit *ytFit = new LeastSquaresFit(fNumPointsSel, tprosel2, yprosel2, erryprosel2);
        vY = ytFit->GetSlope();
        wY = ytFit->GetOffset();
        delete ytFit;
    }
    
    if ( useMFplane ) {
        MedianFit *xtFit = new MedianFit(fNumHitsSel, tprosel3, xprosel3);
        vX = xtFit->GetSlope();
        wX = xtFit->GetOffset();
        delete xtFit;
        MedianFit *ytFit = new MedianFit(fNumHitsSel, tprosel3, yprosel3);
        vY = ytFit->GetSlope();
        wY = ytFit->GetOffset();
        delete ytFit;
    }

    if ( useHFplane ) {
        Float_t ex = 0.001;
        Float_t ey = 0.001;
        Float_t et = 2.5*0.01;
        HoughFit *xhf1 = new HoughFit( tpro, xpro, cpro, 0,et,ex);
    	vX = xhf1->GetSlope();
        wX = xhf1->GetOffset();
    	delete xhf1;
	HoughFit *yhf1 = new HoughFit( tpro, ypro, cpro, 0,et,ey);
    	vY = yhf1->GetSlope();
    	wY = yhf1->GetOffset();
	delete yhf1;
    }

    // calculate normal versor to the TDP
    fNorm.SetXYZ(vY,-vX,wY*vX-wX*vY); 
    fNorm.SetMag(1.);	
    if ( fNorm.Z() > 0 ) fNorm *= -1;

    //use the following only to debug (to see how the projection on plane Z=1 versus time appear after selection) 
    if( fDoGraphUseShape ) {
        Double_t xprog2[fNumPointsSel],yprog2[fNumPointsSel],tprog2[fNumPointsSel];
        Double_t xline2[fNumPointsSel],yline2[fNumPointsSel];
	//Double_t xline2up[fNumPoints],yline2up[fNumPoints];
	//Double_t xline2down[fNumPoints],yline2down[fNumPoints];
	for( Int_t i=0; i<fNumPointsSel; i++ ) {
	    TVector3 dummy = fData.fSpPosSel[i];
	    //xline2[i] = vX * fData.fTimeSel[i] + wX;
	    //yline2[i] = vY * fData.fTimeSel[i] + wY;
	    xline2[i] = vX * tprosel2[i] + wX;
	    yline2[i] = vY * tprosel2[i] + wY;
	    //xline2up[i] = xline2[i] + fMulti2 * a2X;
	    //xline2down[i] = xline2[i] - fMulti2 * a2X;
	    //yline2up[i] = yline2[i] + fMulti2 * a2Y;
	    //yline2down[i] = yline2[i] - fMulti2 * a2Y;
	    xprog2[i] = xprosel2[i]; // xprog2[i] = dummy.X()/dummy.Z();
	    yprog2[i] = yprosel2[i]; // yprog2[i] = dummy.Y()/dummy.Z();
	    tprog2[i] = tprosel2[i]; // tprog2[i] = fData.fTimeSel[i];
	}
	TGraph *gxypro2 = new TGraph(fNumPointsSel,yprog2,xprog2);
	gxypro2->SetNameTitle("gxypro2","xprog2 vs yprog2");
	gxypro2->SetMarkerSize(.5);
	gxypro2->SetMarkerStyle(23);
	TGraph *gxpro2 = new TGraph(fNumPointsSel,tprog2,xprog2);
	gxpro2->SetNameTitle("gxpro2","xprog2 vs time2");
	gxpro2->SetMarkerSize(.5);
	gxpro2->SetMarkerStyle(23);
	TGraph *gypro2 = new TGraph(fNumPointsSel,tprog2,yprog2);
	gypro2->SetNameTitle("gypro2","yprog2 vs time2");
	gypro2->SetMarkerSize(.5);
	gypro2->SetMarkerStyle(23);
	TGraph *gxproretta2 = new TGraph(fNumPointsSel,tprog2,xline2);
	gxproretta2->SetMarkerColor(4);gxproretta2->SetLineColor(4);
	TGraph *gyproretta2 = new TGraph(fNumPointsSel,tprog2,yline2);
	gyproretta2->SetMarkerColor(4);gyproretta2->SetLineColor(4);
	/*TGraph *gxproretta2up = new TGraph(fNumPointsSel,tprog2,xline2up);
	gxproretta2up->SetMarkerColor(2);gxproretta2up->SetLineColor(2);
	TGraph *gyproretta2up = new TGraph(fNumPointsSel,tprog2,yline2up);
	gyproretta2up->SetMarkerColor(2);gyproretta2up->SetLineColor(2);
	TGraph *gxproretta2down = new TGraph(fNumPointsSel,tprog2,xline2down);
	gxproretta2down->SetMarkerColor(2);gxproretta2down->SetLineColor(2);
	TGraph *gyproretta2down = new TGraph(fNumPointsSel,tprog2,yline2down);
	gyproretta2down->SetMarkerColor(2);gyproretta2down->SetLineColor(2);*/
	TCanvas *cUSFP2=new TCanvas("cUSFP2","cUSFP2",1);
	cUSFP2->Divide(3,1);
	cUSFP2->cd(1);gxpro2->Draw("AP");gxproretta2->Draw("PLsame");//gxproretta2up->Draw("PLsame");gxproretta2down->Draw("PLsame");
	cUSFP2->cd(2);gypro2->Draw("AP");gyproretta2->Draw("PLsame");//gyproretta2up->Draw("PLsame");gyproretta2down->Draw("PLsame");
	cUSFP2->cd(3);gxypro2->Draw("AP");
	cUSFP2->Write();
    }
    
    // clear of containers
    cpro.clear(); xpro.clear(); ypro.clear(); tpro.clear();
    errxpro.clear(); errypro.clear();
    xprosel.clear(); yprosel.clear(); tprosel.clear();
    errxprosel.clear(); erryprosel.clear();
    cprosel2.clear(); xprosel2.clear(); yprosel2.clear();
    errxprosel2.clear(); erryprosel2.clear(); tprosel2.clear();
    
    Msg(EsafMsg::Info) << "UseShapeandFindPlane(): " << fNumPointsSel << " pixel and " << fNumHitsSel << " hits selected " << MsgDispatch;

    return;
}

//______________________________________________________________________________
 void TrackDirection2Module::AA1() {
    //
    // ANALITIC APPROXIMATED 1 method.
    // Shower constant angular velocity of the shower approximation.
    // This method initialize the parameters for other all fit methods.
    //

    Msg(EsafMsg::Info) << "Using method AA1()" << MsgDispatch;
	
    fData.fMedDir.SetXYZ(0,0,0);
    fData.fMedTime = 0;
    fData.fMedAlpha = 0;
    vector<Double_t> erralpha;
    vector<Int_t> counts;
    //vectors of data for median fit
    vector<Double_t> alphahits;
    vector<Double_t> timehits;
    Int_t numhits = 0;
	
    for( Int_t i=0; i<fNumPoints; i++ ) {
    	TVector3 dummy = fData.fSpPos[i];
	Double_t alphatmp = ACos(dummy*fW);
	fData.fAlpha.push_back(alphatmp);
	erralpha.push_back(fErrAngle/Sqrt((Double_t)fData.fHits[i]));
	counts.push_back(fData.fHits[i]);
        numhits += fData.fHits[i];
	fData.fMedTime += fData.fTime[i]*fData.fHits[i];
	fData.fMedDir += fData.fSpPos[i]*fData.fHits[i];
        //add hits for median fit
        for( Int_t j(0); j<fData.fHits[i]; j++ ) {
            alphahits.push_back(alphatmp);
            timehits.push_back(fData.fTime[i]);
        }
    }

    fData.fApparentTimeLenght = (Int_t)((fData.fTime[fNumPoints-1]-fData.fTime[0])/fGtuLength);
    fData.fMedTime = fData.fMedTime/fNumHits;
    fData.fMedDir = (fData.fMedDir*(1/(Double_t) fNumHits)).Unit();
    fData.fMedAlpha = ACos(fData.fMedDir*fW);
	
    // shower angular velocity fit: least squares, median and hough 
    Double_t fAngularSpeedLF(0), fQAngularSpeedLF(0), fAngularSpeedtodrawLF(0);
    Double_t fAngularSpeedMF(0), fQAngularSpeedMF(0), fAngularSpeedtodrawMF(0);
    Double_t fAngularSpeedHF(0), fQAngularSpeedHF(0), fAngularSpeedtodrawHF(0);

    if ( fDebugInfo || fDoGraphUseShape || useLFaa1 ) {
        LeastSquaresFit *ASl = new LeastSquaresFit(fNumPoints, fData.fTime, fData.fAlpha, erralpha);
        fAngularSpeedLF  = ASl->GetSlope()*0.01; // uom fix with C()
        fQAngularSpeedLF = ASl->GetOffset();
        fAngularSpeedtodrawLF = ASl->GetSlope();
        delete ASl;
    }
	
    if ( fDebugInfo || fDoGraphUseShape || useMFaa1 ) {
        MedianFit *ASm = new MedianFit(numhits, timehits, alphahits); //CHANGED
        fAngularSpeedMF  = ASm->GetSlope()*0.01; // uom fix with C()
        fQAngularSpeedMF = ASm->GetOffset();
        fAngularSpeedtodrawMF = ASm->GetSlope();
        delete ASm;
    }
		
    if ( fDebugInfo || fDoGraphUseShape || useHFaa1 ) {
        Float_t et = 0.01;       // errors for the hough fit
        Float_t ea = fErrAngle;
        HoughFit *ASh = new HoughFit( fData.fTime, fData.fAlpha, counts, 0,et,ea);
        fAngularSpeedHF = ASh->GetSlope()*0.01; // uom fix with C()
        fQAngularSpeedHF = ASh->GetOffset();
        fAngularSpeedtodrawHF = ASh->GetSlope();
        delete ASh;
    }

    if (fDebugInfo)
        Msg(EsafMsg::Info) << "method AA1(): fAngularSpeedLF = " << fAngularSpeedLF << " fAngularSpeedMF = " << fAngularSpeedMF << " fAngularSpeedHF = " << fAngularSpeedHF << MsgDispatch;

    // use the following only to debug (plot angles alphai vs time) 
    Double_t alphai[fNumPoints],timei[fNumPoints],alphailinelf[fNumPoints],alphailinemf[fNumPoints],alphailinehf[fNumPoints];
    if( fDoGraphUseShape ){
        for( Int_t i=0; i<fNumPoints; i++ ) {
	    alphai[i]=fData.fAlpha[i]*DEG;
	    timei[i]=fData.fTime[i];
	    alphailinelf[i]=(fData.fTime[i]*fAngularSpeedtodrawLF+fQAngularSpeedLF)*DEG;
	    alphailinemf[i]=(fData.fTime[i]*fAngularSpeedtodrawMF+fQAngularSpeedMF)*DEG;
	    alphailinehf[i]=(fData.fTime[i]*fAngularSpeedtodrawHF+fQAngularSpeedHF)*DEG;
	}	
        TGraph *gang = new TGraph(fNumPoints,timei,alphai);
        gang->SetNameTitle("gang","alphai vs time;time;alpha");
        gang->SetMarkerSize(.5);gang->SetMarkerStyle(23);
        TGraph *ganglinelf = new TGraph(fNumPoints,timei,alphailinelf);ganglinelf->SetLineColor(4);
        TGraph *ganglinemf = new TGraph(fNumPoints,timei,alphailinemf);ganglinemf->SetLineColor(6);
        TGraph *ganglinehf = new TGraph(fNumPoints,timei,alphailinehf);ganglinehf->SetLineColor(8);
        TCanvas *cAA1=new TCanvas("cAA1","cAA1",1);cAA1->Divide(1,1);
        cAA1->cd(1);
        gang->Draw("AP");
        ganglinelf->Draw("PLsame");
        ganglinemf->Draw("PLsame");
        ganglinehf->Draw("PLsame"); 
        TLegend *legend = new TLegend(0.8,0.9,1,0.6);
        legend->AddEntry(ganglinelf,"least squares fit","l");
        legend->AddEntry(ganglinemf,"median fit","l");
        legend->AddEntry(ganglinehf,"hough fit","l");
        legend->Draw();
        cAA1->Write();
        delete cAA1;
        delete gang;
        delete ganglinelf;
        delete ganglinemf;
        delete ganglinehf;
        delete legend;
    }

    // selection of fit method
    if ( useLFaa1 ) {
        fAngularSpeed = fAngularSpeedLF;
    } else if ( useMFaa1 ) {
        fAngularSpeed = fAngularSpeedMF;
    } else if ( useHFaa1 ) {
        fAngularSpeed = fAngularSpeedHF;
    } else {
        fAngularSpeed = fAngularSpeedLF;
        Msg(EsafMsg::Warning) << "AA1() fith method non specified in config file. Saving least squares fit results." << MsgDispatch;
    }

    // check if the result of the fit is a 'not a number'
    if ( IsNaN(fAngularSpeed) ) {
        Msg(EsafMsg::Warning) << "NaN angular speed, skipping this event" << MsgDispatch;
        fAA1done = kTRUE;
        fAA1nan = kTRUE;
        return;
    }
	
    fHmax = 5; // HMax initialized to 5 kilometers
    fRmax = CalculateRmax(fHmax);

    if ( fDebugInfo ) { // if debug display results for all fith methods
        fBetaInit = 2*ATan(CLUCE/(fAngularSpeedLF*fRmax)) - fData.fMedAlpha;
        fBeta = PI - fBetaInit;		
        while(fBeta > 2*PI) fBeta = fBeta - 2*PI;
        CalculatefromBetaEASdir();
        Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirLF = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch; 
	
        fBetaInit = 2*ATan(CLUCE/(fAngularSpeedMF*fRmax)) - fData.fMedAlpha;
        fBeta = PI - fBetaInit;
        while(fBeta> 2*PI) fBeta = fBeta - 2*PI;
        CalculatefromBetaEASdir();
        Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirMF = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch; 
	
        fBetaInit = 2*ATan(CLUCE/(fAngularSpeedHF*fRmax)) - fData.fMedAlpha;
        fBeta = PI - fBetaInit;		
        while(fBeta > 2*PI)	fBeta = fBeta - 2*PI;
        CalculatefromBetaEASdir();
        Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirHF = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch;
    }
    
    // first value of beta
    fBetaInit = 2*ATan(CLUCE/(fAngularSpeed*fRmax)) - fData.fMedAlpha;
    fBeta = PI - fBetaInit;	
    while(fBeta > 2*PI)	fBeta = fBeta - 2*PI;
    CalculatefromBetaEASdir();

    // recalculate Hmax
    fHmax = FindHmax( fTHETAreco );
    fRmax = CalculateRmax( fHmax );
    TVector3 nmax = fData.fMedDir;
    fData.VectorRmax.SetMagThetaPhi(fRmax,nmax.Theta(),nmax.Phi());
	
    // recalculate beta
    fBetaInit = 2*atan(CLUCE/(fAngularSpeed*fRmax)) - fData.fMedAlpha;
    fBeta = PI - fBetaInit;
    while(fBeta > 2*PI)	fBeta = fBeta - 2*PI;
    CalculatefromBetaEASdir();
	
    fHmax = FindHmax(fTHETAreco);
    Msg(EsafMsg::Info) << "method AA1(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch; 
	
    fAA1done = kTRUE;
}

//______________________________________________________________________________
 void TrackDirection2Module::AA2() {
    //
    // ANALYTIC APPROXIMATED 2 method.
    // Approximation: Shower velocity on a plane perpendicular to the detector
    //                axis is constant
    //
	
    Msg(EsafMsg::Info) << "Using method AA2()" << MsgDispatch;
    if ( !fAA1done ) AA1();

    if ( fAA1nan ) return;
    
    vector<Double_t> xprov,yprov,errxyprov;
    vector<Double_t> errxprov,erryprov;
    for(Int_t i=0;i<fNumPoints;i++){
	TVector3 dummy = fData.fSpPos[i];
	xprov.push_back(dummy.x()*(fHISS-fHmax)/dummy.z());
	yprov.push_back(dummy.y()*(fHISS-fHmax)/dummy.z());
        Double_t dx = (fHISS-fHmax) * DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
        Double_t dy = (fHISS-fHmax) * DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
	errxyprov.push_back(1./Sqrt((Double_t)fData.fHits[i]));
	errxprov.push_back( dx / fData.fHits[i]);
	erryprov.push_back( dy / fData.fHits[i]);
    }
	
    Double_t speedx,speedy;
    LeastSquaresFit *fitvx = new LeastSquaresFit(fNumPoints, fData.fTime, xprov, errxprov);
    speedx = fitvx->GetSlope()*0.01; // why?
    delete fitvx;
    LeastSquaresFit *fitvy = new LeastSquaresFit(fNumPoints, fData.fTime, yprov, erryprov);
    speedy = fitvy->GetSlope()*0.01; // why?
    delete fitvy;
	
    Double_t speed = Sqrt(speedx*speedx + speedy*speedy);
    Double_t gthetamax = fData.fMedAlpha;
    Double_t betacontrol = fBeta;
    Double_t gammaV = 2*ATan((speed/CLUCE) * Sin(gthetamax));
    Double_t betaV;
	
    if(betacontrol*DEG>(90-gthetamax*DEG) && betacontrol*DEG<90){
	betaV = -gammaV + gthetamax;
    } else if (betacontrol*DEG > (90-gthetamax*DEG) && betacontrol*DEG > 90 ) {
	    betaV = gthetamax + gammaV;		
    }
	
    fBeta = betaV;
    while(fBeta > 2*PI) fBeta=fBeta-2*PI;
    CalculatefromBetaEASdir();	
    fHmax = FindHmax(fTHETAreco);
    Msg(EsafMsg::Info) << "method AA2(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch; 	
}

//______________________________________________________________________________
 void TrackDirection2Module::NE1() {
    //
    // NUMERICAL EXACT 1 method
    // Chi-square minimization of the difference between arrival times of photons
    // measured and teoretically computed
    //
    
    Msg(EsafMsg::Info) << "using method NE1()" << MsgDispatch;
    if( !fAA1done ) AA1();
    if ( fAA1nan ) return;

    TMinuit *minuitNE1 = new TMinuit(3);	
    minuitNE1->SetObjectFit(&fData);	
    minuitNE1->SetPrintLevel(fMinuitOutputLevel);

    TString  namebeta="beta", nametmax="tmax", namehmax="hmax";
    Int_t ierflg(0);
    Double_t fcnout,edm,errdef;
    Int_t nvpar,nparx,icstat;
    Double_t value, errv, vlow,vup;
    Int_t gmaxcall=1000;
    Double_t deltabeta,deltatmax;
    Double_t betaStart = PI - fBeta;
    Double_t hmaxStart = fHmax;
    Double_t tmaxStart = fData.fMedTime*100;
    Double_t step[3] = {0.2,2.500,0.3}; //about 10 deg (beta), 2500ns (tmax), 0.3 km (hmax)
    Double_t arglist[10];
	
    minuitNE1->SetFCN(ChiSquareNE1);
    minuitNE1->mnparm(0, "beta",betaStart , step[0], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER beta.");
    minuitNE1->mnparm(1, "tmax",tmaxStart , step[1], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER tmax.");
    minuitNE1->mnparm(2, "hmax",hmaxStart , step[2], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER hmax.");
	
    if( fFixTmaxNumeric ) {
	arglist[0] = 2; 
	minuitNE1->mnexcm("FIX",arglist,1,ierflg); if (ierflg) Printf(" .....UNABLE to fix parameter tmax");	
    }
    arglist[0] = 3; 
    minuitNE1->mnexcm("FIX",arglist,1,ierflg); if (ierflg) Printf(" .....UNABLE to fix parameter tmax");	
    arglist[0] = gmaxcall;
    minuitNE1->mnexcm("MIGRAD",arglist ,1 , ierflg);if (ierflg) 	Printf(" ......UNABLE to minimize with MIGRAD");
	
    minuitNE1->mnstat(fcnout,edm,errdef,nvpar,nparx,icstat);
    minuitNE1->mnpout(0,namebeta, value,errv,vlow,vup,ierflg);
    fBeta = PI - value;
    deltabeta = errv;
	
    minuitNE1->mnpout(2,nametmax, value,errv,vlow,vup,ierflg);
    fTmaxFit = value;
    deltatmax = errv;
    CalculatefromBetaEASdir();		
    fHmax=FindHmax(fTHETAreco);
    Msg(EsafMsg::Info) << "method NE1(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch; 
}

//______________________________________________________________________________
void ChiSquareNE1(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag) { 
    //
    // Function FCN called by minuit for the algorithm NE1
    //
	
    chisqnorm = 0;	
    ContainerData *fDataMin = (ContainerData*) gMinuit->GetObjectFit();
	
    Int_t numhits(0);
    Double_t chisq(0);
    Double_t gerrort = fDataMin->fGtuLength;
    Double_t Rm = (fDataMin->VectorRmax).Mag();
    Double_t Am = fDataMin->fMedAlpha;
    for(Int_t i=0;i<fDataMin->fNumPoints;i++){  
        Double_t TimeExpected = ( par[1] - Rm * ( ( Sin(Am-fDataMin->fAlpha[i]) + Sin(fDataMin->fAlpha[i]+par[0]) - Sin(Am+par[0]) )/ Sin(fDataMin->fAlpha[i]+par[0]) )/CLUCE);
	chisq += Power((fDataMin->fTime[i]*100-TimeExpected)/gerrort,2)*fDataMin->fHits[i]; 
	numhits += fDataMin->fHits[i];
    }
	
    chisqnorm = (Double_t)chisq/numhits;
}


//______________________________________________________________________________
 void TrackDirection2Module::NE2() {
    //
    // NUMERICAL EXACT 2 method
    // Chi-square minimization of angle between the versors of pixels in FOV
    // and the corresponding vectors from points of the track and the detector
    //

    Msg(EsafMsg::Info) << "using method NE2()" << MsgDispatch;
    if( !fAA1done ) AA1();
    if ( fAA1nan ) return;

    TMinuit *minuitNE2 = new TMinuit(3);	
    minuitNE2->SetObjectFit(&fData);   
    minuitNE2->SetPrintLevel(fMinuitOutputLevel);

    TString nametheta="theta",namephi="phi",nametmax="tmax";
    Int_t ierflg(0);
    Double_t fcnout,edm,errdef;
    Int_t nvpar,nparx,icstat;
    Double_t value, errv, vlow,vup;
    Int_t gmaxcall = 1000;
    Double_t deltatheta,deltaphi,deltatmax;
	
    Double_t thetastart = fTHETAreco;
    Double_t phistart = fPHIreco;
    Double_t tmaxStart = fData.fMedTime*100;  //microseconds
    Double_t step[3] = {0.2,0.2,2.500};
    Double_t arglist[10];
	
    minuitNE2->SetFCN(ChiSquareNE2);
    minuitNE2->mnparm(0, "theta",thetastart , step[0], 0,0,ierflg); if (ierflg)	Printf("......UNABLE TO DEFINE PARAMETER theta.");
    minuitNE2->mnparm(1, "phi",phistart , step[1], 0,0,ierflg); if (ierflg) Printf(".....UNABLE TO DEFINE PARAMETER phi.");
    minuitNE2->mnparm(2, "tmax",tmaxStart , step[2], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER tmax.");
	
    if( fFixTmaxNumeric ) {
	arglist[0] = 3; 
	minuitNE2->mnexcm("FIX",arglist,1,ierflg);if (ierflg) 	Printf(" .....UNABLE to fix parameter tmax");	
    }
    arglist[0] = gmaxcall;		
    minuitNE2->mnexcm("MIGRAD",arglist ,1 , ierflg);if (ierflg) 	Printf(" .......UNABLE to minimize with migrad");
				
    minuitNE2->mnstat(fcnout,edm,errdef,nvpar,nparx,icstat);
    minuitNE2->mnpout(0,nametheta, value,errv,vlow,vup,ierflg);
    fTHETAreco = value;
    deltatheta = errv;
		
    minuitNE2->mnpout(1,namephi, value,errv,vlow,vup,ierflg);
    fPHIreco = value;
    deltaphi = errv;
		
    minuitNE2->mnpout(2,nametmax, value,errv,vlow,vup,ierflg);
    fTmaxFit = value;
    deltatmax = errv;
	
    fEASDir.SetMagThetaPhi(1,fTHETAreco,fPHIreco);
    fTHETAreco = fEASDir.Theta();
    fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*PI) );
    fDeltaTheta = fTHETAreco-fTrueTheta;
    fDeltaPhi = fPHIreco-fTruePhi;
    fDeltaEASDir = fEASDir.Angle(fTrueDir);
	
    fHmax = FindHmax(fTHETAreco);
    Msg(EsafMsg::Info) << "method NE2(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch; 
}

//______________________________________________________________________
void ChiSquareNE2(Int_t &npar, Double_t *gin, Double_t &chisqnorm, Double_t *par, Int_t iflag) {
    //
    // Function FCN called by minuit for the algorithm NE2
    //
	
    ContainerData *fDataMin = (ContainerData*) gMinuit->GetObjectFit();
	
    Int_t numhits(0);
    Double_t PidotRi(0), modRi(0), Li(0), betaprimo(0);
    TVector3 dummyRmax = fDataMin->VectorRmax;
    Double_t Xm = dummyRmax.x();
    Double_t Ym = dummyRmax.y();
    Double_t Zm = dummyRmax.z();
    Double_t Rm = dummyRmax.Mag();
    Double_t gerrorpix = 0.0018; //rad, about 0.1 deg
    Double_t sctmp = Sin(par[0])*Cos(par[1]-PI);
    Double_t sstmp = Sin(par[0])*Sin(par[1]-PI);
    Double_t ctmp = -Cos(par[0]);
    Double_t thetaprimotmp = ACos((-Xm*sctmp-Ym*sstmp-Zm*ctmp)/Rm);
	
    chisqnorm=0;
    for(Int_t i=0; i<fDataMin->fNumPoints; i++){  
        betaprimo = 2*ATan(1/((1/Tan(thetaprimotmp/2))+CLUCE*(fDataMin->fTime[i]*100-par[2])/(Rm*Sin(thetaprimotmp))));
	Li = Rm*(Cos(thetaprimotmp) - Sin(thetaprimotmp) / Tan(betaprimo));
	PidotRi = fDataMin->fSpPos[i].x()*(Xm+Li*sctmp) + fDataMin->fSpPos[i].y()*(Ym+Li*sstmp) + fDataMin->fSpPos[i].z()*(Zm+Li*ctmp);
	modRi = Power(Xm+Li*sctmp,2) + Power(Ym+Li*sstmp,2) + Power(Zm+Li*ctmp,2);
	chisqnorm += Power((ACos(PidotRi/Sqrt(modRi)))/gerrorpix,2)*fDataMin->fHits[i]; 
	numhits += fDataMin->fHits[i];
    }	
    chisqnorm = (Double_t)chisqnorm/numhits;
}

//______________________________________________________________________________
 void TrackDirection2Module::AE1() {	
    //
    // ANALYTICAL EXACT 1 method
    // Fit using exact relations between pixel directions in FOV and photons
    // arrival times. This method doesn't require the knowledge of the TDP.
    //
    
    Msg(EsafMsg::Info) << "Using method AE1()" << MsgDispatch;
    if( !fAA1done ) AA1();
    if ( fAA1nan ) return;
	
    vector<Double_t> alpha;
    vector<Double_t> Ri;
    vector<Double_t> dLi;
    vector<Double_t> timestart;
    vector<Double_t> nh;
    vector<Double_t> x;
    vector<Double_t> y;
    vector<Double_t> z;
    vector<Double_t> ex;
    vector<Double_t> ey;
    vector<Double_t> ez;
    Double_t gnthetamax = fData.fMedDir.Theta();
    Double_t gnphimax = fData.fMedDir.Phi();
    Double_t timemax = fData.fMedTime*100; // microseconds
	
    Int_t pointsAE1(0),hitsAE1(0);
    Double_t Ritmp(0), dLitmp(0);
    for ( Int_t i=0 ; i<fNumPoints; i++ ){
    	TVector3 dummy = fData.fSpPos[i];
	Double_t dummyTime = fData.fTime[i]*100;//mus
        Double_t cost = Cos(dummy.Theta())*Cos(gnthetamax) + Sin(dummy.Theta())*Sin(gnthetamax)*Cos(dummy.Phi()-gnphimax);
	if( Abs(dummyTime-timemax)/fGtuLength>1 && Abs(dummyTime-timemax)/fGtuLength>fData.fApparentTimeLenght/10 && fData.fHits[i]>2 ){
	    if( dummyTime > timemax){
	        Ritmp = (2*CLUCE*(dummyTime-timemax)*fRmax + Power(CLUCE*(dummyTime-timemax),2))/(2*(fRmax*(1-cost) 
                    + CLUCE*(dummyTime-timemax)));
		dLitmp = -Sqrt( Ritmp*Ritmp + fRmax*fRmax - 2*Ritmp*fRmax*cost );
	    }
	    if(dummyTime < timemax){
		Ritmp = (2*CLUCE*(-dummyTime+timemax)*fRmax - Power(CLUCE*(dummyTime-timemax),2))/(2*(fRmax*(-1+cost) 
                    + CLUCE*(-dummyTime+timemax)));
		dLitmp = Sqrt( Ritmp*Ritmp + fRmax*fRmax - 2*Ritmp*fRmax*cost );
	    }	
	    alpha.push_back(fData.fAlpha[i]);
	    Ri.push_back(Ritmp);
	    dLi.push_back(dLitmp);
	    timestart.push_back(-dLitmp/CLUCE);
	    nh.push_back(fData.fHits[i]);
	    x.push_back(Ritmp*Sin(dummy.Theta())*Cos(dummy.Phi()));
	    ex.push_back(0.1/fData.fHits[i]);
            y.push_back(Ritmp*Sin(dummy.Theta())*Sin(dummy.Phi()));
            ey.push_back(0.1/fData.fHits[i]);
	    z.push_back(Ritmp*Cos(dummy.Theta()));
	    ez.push_back(0.1/fData.fHits[i]);
	    hitsAE1 += fData.fHits[i];
	    pointsAE1++;
	}
    }
	
    alpha.push_back(fData.fMedAlpha);
    dLi.push_back(0);
    timestart.push_back(0);
    Ri.push_back(fRmax);
    x.push_back(fData.VectorRmax.x());
    y.push_back(fData.VectorRmax.y());
    z.push_back(fData.VectorRmax.z());
	
    Double_t vx,vy,vz;
    LeastSquaresFit *fitvx = new LeastSquaresFit(pointsAE1, timestart, x, ex);
    vx = fitvx->GetSlope();
    delete fitvx;
    LeastSquaresFit *fitvz = new LeastSquaresFit(pointsAE1, timestart, z, ez);
    vz = fitvz->GetSlope();	
    delete fitvz;
    LeastSquaresFit *fitvy = new LeastSquaresFit(pointsAE1, timestart, y, ey);
    vy = fitvy->GetSlope();	
    delete fitvy;

    if ( IsNaN(vx) ||IsNaN(vy) || IsNaN(vz) ) {
        Msg(EsafMsg::Warning) << "AE1() fit returns a NaN, skipping" << MsgDispatch;
        return;
    }
    
    Double_t m1 = vx/vz;
    Double_t m2 = vy/vz;
    Double_t tT = Sqrt(m1*m1+m2*m2);
    Double_t tP = m2/m1;
	
    fTHETAreco = ATan(tT);
    fPHIreco = ATan(tP);
    if( tP>0 ){
	fPHIreco = ( m2>0 ? atan(tP) : PI+atan(tP) );
    } else {
	fPHIreco = ( m2>0 ? PI+atan(tP) : (2*PI)+atan(tP) );
    }
	
    fEASDir.SetMagThetaPhi(1,fTHETAreco,fPHIreco);
    fTHETAreco = fEASDir.Theta();
    fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*PI) );
    fDeltaTheta = fTHETAreco - fTrueTheta;
    fDeltaPhi = fPHIreco - fTruePhi;
    fDeltaEASDir = fEASDir.Angle(fTrueDir);	
	
    fHmax = FindHmax(fTHETAreco);
    Msg(EsafMsg::Info) << "method AE1(): fDeltaEASDir = " <<fDeltaEASDir*DEG << " deg ( dT = " << fDeltaTheta*DEG << " , dP = " << fDeltaPhi*DEG << " )" << MsgDispatch; 

}
//______________________________________________________________________________
 void TrackDirection2Module::all() {	
    //
    // Execute all methods for reconstructing the shower direction.
    //
    
    Msg(EsafMsg::Info) << "Using all methods .. " << MsgDispatch;

    AA1();
    fTHETArecoAA1 = fTHETAreco;
    fPHIrecoAA1 = fPHIreco;
    fDeltaEASDirAA1 = fDeltaEASDir;
    AA2();
    fTHETArecoAA2 = fTHETAreco;
    fPHIrecoAA2 = fPHIreco;
    fDeltaEASDirAA2 = fDeltaEASDir;
    NE1();
    fTHETArecoNE1 = fTHETAreco;
    fPHIrecoNE1 = fPHIreco;
    fDeltaEASDirNE1 = fDeltaEASDir;
    NE2();
    fTHETArecoNE2 = fTHETAreco;
    fPHIrecoNE2 = fPHIreco;
    fDeltaEASDirNE2 = fDeltaEASDir;
    AE1();
    fTHETArecoAE1 = fTHETAreco;
    fPHIrecoAE1 = fPHIreco;
    fDeltaEASDirAE1 = fDeltaEASDir;
	
    //if all is chosen are saved only AA1 results
    fTHETAreco = fTHETArecoAA1;
    fPHIreco = fPHIrecoAA1;
}

//______________________________________________________________________________
 void TrackDirection2Module::CalculatefromBetaEASdir() {
    //
    // Calculate vector of EAS direction using angle fBeta 
    // and the equation of TrackDirectionPlane
    //

    fEASDir = Cos(fBeta) * fW + Sin(fBeta) * fU;
    fTHETAreco = fEASDir.Theta();
    fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*PI) );
    fDeltaTheta = fTHETAreco - fTrueTheta;
    fDeltaPhi = fPHIreco - fTruePhi;
    fDeltaEASDir = fEASDir.Angle(fTrueDir);
    return; 
}

//______________________________________________________________________________
 Double_t TrackDirection2Module::FindHmax( Double_t theta ) {
    //
    // Find Hmax with the Linsley parametrization of the atmosphere 
    // depending on the zenith angle of the shower and the value of Xmax
    // (fixed value is only a first approximation).
    //

    Double_t Xmax = 831; //g/cm^2
    Double_t thetaloc = theta; //it should be zenith angle in local reference frame, to improve.
    Double_t Xv = Xmax*Cos(thetaloc);
	
    Double_t X0 = 1036.1; // g/cm^2
    Double_t X4 = 631.1;  // g/cm^2
    Double_t X10 = 271.1; // g/cm^2
    Double_t al(0), bl(0), cl(0);
    Double_t val(0);
    if( Xv<X0 && Xv>X4 ){
        cl = 9.9418638;  // km
	bl = 1222.6562;  // g/cm2
	al = -186.5562;  // g/cm2
    } else if( Xv<X4 && Xv>X10 ) {
	cl = 8.7815355;  // km
	bl = 1144.9069;  // g/cm2
	al = -94.9199;   // g/cm2
    } else if( Xv<X10 ) {
	cl = 6.3614304;  // km
	bl = 1305.5948;  // g/cm2
	al = 0.61289;    // g/cm2
    } else {
	Msg(EsafMsg::Warning) <<"...something was wrong in TrackDirection2Module::FindHmax"<< MsgDispatch;
	return 0;
    }
    val = -cl * Log((Xmax * Cos(thetaloc)-al)/bl);
    return val;
}

//______________________________________________________________________________
 Double_t TrackDirection2Module::CalculateRmax( Double_t hmax ) {
    //
    // Method to geometrically find Rmax using Hmax and the versor pointing 
    // to the maximum of the shower
    //
    
    Double_t thmax = fData.fMedDir.Theta();
    Double_t rmax = (RT+fHISS)*TMath::Cos(thmax) - TMath::Sqrt(TMath::Power(RT+hmax,2) - TMath::Power((RT+fHISS)*TMath::Sin(thmax),2));
    return rmax;
}

//______________________________________________________________________________
 Double_t TrackDirection2Module::DeltaX( Double_t theta, Double_t phi, Double_t errtheta, Double_t errphi ) {
    //
    // Calculate error on X projection on the focal surface of a given point
    // on the unitary sphere
    //

    return Sqrt( Power((Cos(phi)*Cos(theta)+Sin(theta))*errtheta/(Cos(theta)*Cos(theta)),2.) +
                 Power(Sin(theta)*Sin(phi)*errphi/Cos(theta),2.));
}

//______________________________________________________________________________
 Double_t TrackDirection2Module::DeltaY( Double_t theta, Double_t phi, Double_t errtheta, Double_t errphi ) {
    //
    // Calculate error on Y projection on the focal surface of a given point
    // on the unitary sphere
    //

    return Sqrt( Power((Sin(phi)*Cos(theta)+Sin(theta))*errtheta/(Cos(theta)*Cos(theta)),2.) +
                Power(Sin(theta)*Cos(phi)*errphi/Cos(theta),2.));
}

//___________________________________________________
void ContainerData::Clear() {
    //
    // Clear of the container data class
    //
    
    fPixXYZ.clear();
    fPos.clear();
    fErr.clear();
    fSpPos.clear();
    fSpErr.clear();
    fSpPosSel.clear();
    fTime.clear();
    fTimeSel.clear();
    fAlpha.clear();
    fSigmaPhi.clear();
    fSigmaTheta.clear();
    fSigmaPhiSel.clear();
    fSigmaThetaSel.clear();
    fHits.clear();
    fHitsSel.clear();
    fNumPoints = 0;
    fNumHits = 0;
    fCentroid.SetXYZ(0,0,0);
    fMedDir.SetXYZ(0,0,0);
    VectorRmax.SetXYZ(0,0,0);
    fMedTime = 0;
    fMedAlpha = 0;
    fApparentTimeLenght = 0;
    fGtuLength = 0;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005