#include "IExample2010aAnalysis1Module.hxx"

// set the environment variable 
// CMTEXTRATAGS to OAANALYSIS_WITHOUT_RECON_CALIB
// to disable this code
#ifndef OAANALYSIS_WITHOUT_RECON_CALIB

#include "ICOMETLog.hxx"

#include <IComboHit.hxx>
#include <IG4HitSegment.hxx>
#include "IG4PrimaryVertex.hxx"
#include <IG4Trajectory.hxx>
#include "IG4Trajectory.hxx"
#include <IG4TrajectoryPoint.hxx>
#include <IG4VHit.hxx>
//#include <IHit.hxx>
#include "IMultiHit.hxx"
#include <ICalibManager.hxx>
#include "IMCHit.hxx"
#include "IRealDatum.hxx"
#include "IReconPID.hxx"
#include <IOARuntimeParameters.hxx>
#include "ICOMETContext.hxx"
#include <IFGDChannelCalibrator.hxx>
#include <ITPCChannelCalibrator.hxx>
#include "ITPCReco.hxx"
#include "CLHEP/Matrix/Matrix.h"
#include <TrackingUtils.hxx>

#include "IRecPackManager.hxx"

ClassImp(COMET::IExample2010aAnalysis1Module);
ClassImp(COMET::IExample2010aAnalysis1Module::ITrackFromTPCPids);
ClassImp(COMET::IExample2010aAnalysis1Module::IFGDHit);
ClassImp(COMET::IExample2010aAnalysis1Module::IP0DHit);
ClassImp(COMET::IExample2010aAnalysis1Module::ISMRDHit);
ClassImp(COMET::IExample2010aAnalysis1Module::ITrueNuVertex);

#define CVSTAG "\
  $Name: v5r19 $"
#define CVSID  "\
  $Id: IExample2010aAnalysis1Module.cxx,v 1.8 2013/02/04 12:52:52 db211 Exp $"

COMET::IExample2010aAnalysis1Module::IExample2010aAnalysis1Module(const char *name,
const char *title):
fTracksFromTPCPids(new TClonesArray("COMET::IExample2010aAnalysis1Module::ITrackFromTPCPids", 10)),
fFGDHits(new TClonesArray("COMET::IExample2010aAnalysis1Module::IFGDHit", 200)),
fP0DHits(new TClonesArray("COMET::IExample2010aAnalysis1Module::IP0DHit", 200)),
fSMRDHits(new TClonesArray("COMET::IExample2010aAnalysis1Module::ISMRDHit", 200)),
fTrueNuVertices(new TClonesArray("COMET::IExample2010aAnalysis1Module::ITrueNuVertex", 10)) {
    SetNameTitle(name, title);
    // Enable this module by default:
      fIsEnabled = kFALSE;
    fDescription = "2010a Analysis Example 1: adapted from Claudio Giganti's analysis code";
    fCVSTagName = CVSTAG;
    fCVSID = CVSID;
}


COMET::IExample2010aAnalysis1Module::~IExample2010aAnalysis1Module() {}

Bool_t COMET::IExample2010aAnalysis1Module::ProcessFirstEvent(COMET::ICOMETEvent& event) {
    fTPCRecon = new ITPCReco();
    COMET::ICalibManager::Get().AddCalibrator(new COMET::ITPCChannelCalibrator());
    COMET::ICalibManager::Get().AddCalibrator(new COMET::IFGDChannelCalibrator());
    return true;
}


bool COMET::IExample2010aAnalysis1Module::IsFullEventWorthSaving(COMET::ICOMETEvent& event) {
    return (fNTracksFromTPCPids ? true : false);
}


void COMET::IExample2010aAnalysis1Module::InitializeModule() {

    fTimeOffset = COMET::IOARuntimeParameters::Get().GetParameterD("ReconTPC.Reco.LikFit.TimeOffset");

}


void COMET::IExample2010aAnalysis1Module::InitializeBranches() {

    fOutputTree->Branch("NTPCHits", &fNTPCHits, "NTPCHits/I", fBufferSize);
    fOutputTree->Branch("NTracksFromTPCPids", &fNTracksFromTPCPids, "NTracksFromTPCPids/I", fBufferSize);
    fOutputTree->Branch("NFGDHits", &fNFGDHits, "NFGDHits/I", fBufferSize);
    fOutputTree->Branch("NP0DHits", &fNP0DHits, "NP0DHits/I", fBufferSize);
    fOutputTree->Branch("NSMRDHits", &fNSMRDHits, "NSMRDHits/I", fBufferSize);
    fOutputTree->Branch("NTrueNuVertices", &fNTrueNuVertices, "NTrueNuVertices/I", fBufferSize);

    fOutputTree->Branch("TracksFromTPCPids", "TClonesArray", &fTracksFromTPCPids, fBufferSize, fSplitLevel);
    fOutputTree->Branch("FGDHits", "TClonesArray", &fFGDHits, fBufferSize, fSplitLevel);
    fOutputTree->Branch("P0DHits", "TClonesArray", &fP0DHits, fBufferSize, fSplitLevel);
    fOutputTree->Branch("SMRDHits", "TClonesArray", &fSMRDHits, fBufferSize, fSplitLevel);
    fOutputTree->Branch("TrueNuVertices", "TClonesArray", &fTrueNuVertices, fBufferSize, fSplitLevel);

}


bool COMET::IExample2010aAnalysis1Module::FillTree(COMET::ICOMETEvent& event) {

    // Initialise
    fNTracksFromTPCPids = 0;
    fNFGDHits = 0;
    fNP0DHits = 0;
    fNSMRDHits = 0;
    fNTrueNuVertices = 0;

    fTracksFromTPCPids->Clear();
    fFGDHits->Clear();
    fP0DHits->Clear();
    fSMRDHits->Clear();
    fTrueNuVertices->Clear();

    /* to fill a new output object, for example a new ITrackFromTPCPids object, do
    ITrackFromTPCPids* trackFromTPCPids = new((*fTracksFromTPCPids)[fNTracksFromTPCPids]) ITrackFromTPCPids;
    // fill the variables
    trackFromTPCPids->npoint = 493;

    fNTracksFromTPCPids++;
    */
    //int evtID = event.GetEventId();
    //int runID = event.GetRunId();

    //int rCode = 0;

    //nevt++;

    /*
    char evnum[100]; sprintf(evnum,"%4.1d", evtID);
    char rcode[100]; sprintf(rcode,"%2.1d", rCode);
    TString evstr = TString("Event No.:") + TString(evnum)
        +  TString("   Reaction code:") + TString(rcode);

    COMETDebug( evstr.Data() );
    COMETDebug( " ------------------------ " );

    int trigbit =(int) event.GetHeader().GetTriggerBits();

    COMETDebug(" Trigger bits "<<std::hex<<"0x"<<trigbit<<std::dec);

    */
    //if( trigbit != 0x100 && trigbit !=0x80 ) return true; ///skip pedestals
    //if( trigbit == 0x4 ) return true; ///skip laser

    //tree.nevent=evtID;
    //tree.run=runID;

    //#ifdef MC
    if(event.GetContext().IsMC()) {
        GetTrueVertex(event);
    }
    //#endif

    // Get the FGD hits.
    //COMETDebug("look at FGD");
    COMET::IHandle<COMET::IHitSelection> fgd(event.GetHitSelection("fgd"));
    if( !fgd ) {
        COMET::ICalibManager::Get().Calibrate("fgd");
        COMETWarn( " NO FGD hits found --> Try to calibrate ..." );
    }
    fgd = event.GetHitSelection("fgd");
    if( !fgd ) {
        COMETWarn( " NO FGD hits found after calibration" );
        //tree.ntrack=0;
        //treefile->Fill();
        return true;
    }
    //COMETDebug("fill fgd hits");
    //COMETDebug("look at P0D");
    COMET::IHandle<COMET::IHitSelection> p0d(event.GetHitSelection("p0d"));

    //COMETDebug("look at SMRD");
    COMET::IHandle<COMET::IHitSelection> smrd(event.GetHitSelection("smrd"));

    if(fgd)
        FillFGDHits(fgd);

    if(p0d)
        FillP0DHits(p0d);

    if(smrd)
        FillSMRDHits(smrd);

    // Get the TPC hits.
    COMET::IHandle<COMET::IHitSelection> tpc(event.GetHitSelection("tpc"));
    if( !tpc ) {
        COMET::ICalibManager::Get().Calibrate("tpc");
        COMETWarn( " NO TPC hits found --> Try to calibrate ..." );
    }
    tpc = event.GetHitSelection("tpc");
    if( !tpc ) {
        COMETWarn( " NO TPC hits found after calibration. Skip event !!!" );
        //tree.ntrack=0;
        //treefile->Fill();
        return true;
    }

    COMETDebug( " Number of TPC Hits " << tpc->size() );

    fNTPCHits = tpc->size();

    COMET::IHandle<COMET::IAlgorithmResult> tpcHits(fTPCRecon->HitPreparation());


    // Comment out as functionality of ITPCLikFit has changed
    // fTPCRecon->GetTPCLikFit()->FitSigma(true);

    COMET::IHandle<COMET::IAlgorithmResult> FitResult(fTPCRecon->Process(*tpcHits));

    if( FitResult ) {

        COMET::IHandle<COMET::IReconObjectContainer> TrackCont(FitResult->GetResultsContainer("TPCPids"));

        if( TrackCont ) {

            //FillFGDHits(fgd);

            //COMETDebug(tree.fgdX[0]<<" "<<tree.fgdY[0]<<" "<<tree.fgdZ[0]);
            int itr=0;

                                    // Loop over the received tracks
            for( COMET::IReconObjectContainer::iterator  tr = TrackCont->begin();
            tr < TrackCont->end(); tr++) {
                // Create TPC track object to fill
                ITrackFromTPCPids* trackFromTPCPids = new((*fTracksFromTPCPids)[fNTracksFromTPCPids]) ITrackFromTPCPids;
                double x,y,z,ty,tx; //rho;
                //double ex,ey,etx,ety,ep,ek;
                double etx,ety;
                //double p,q;

                COMET::IHandle<COMET::IReconPID> trackPiD = *tr;
                COMET::IHandle<COMET::IReconObjectContainer>  trackcont = trackPiD->GetConstituents();

                COMET::IHandle<COMET::IReconTrack> track = *(trackcont->begin());

                COMET::IReconNodeContainer &nodes = track->GetNodes();

                trackFromTPCPids->goodness=track->GetQuality();

                if( (*tr)->Get< COMET::IRealDatum >("tpcPid_Craw") ) {
                    trackFromTPCPids->dedxmeas = (*tr)->Get< COMET::IRealDatum >("tpcPid_Craw")->GetValue();
                    trackFromTPCPids->dedxcorr = (*tr)->Get< COMET::IRealDatum >("tpcPid_Ccorr")->GetValue();
                    trackFromTPCPids->nTrun = (*tr)->Get< COMET::IRealDatum >("tpcPid_NTrun")->GetValue();
                    trackFromTPCPids->SampleLength = (*tr)->Get< COMET::IRealDatum >("tpcPid_SampleLength")->GetValue();

                    trackFromTPCPids->dedxexpele = (*tr)->Get< COMET::IRealDatum >("tpcPid_dEdxexpEle")->GetValue();
                    trackFromTPCPids->dedxexpmuon = (*tr)->Get< COMET::IRealDatum >("tpcPid_dEdxexpMuon")->GetValue();
                    trackFromTPCPids->dedxexppion = (*tr)->Get< COMET::IRealDatum >("tpcPid_dEdxexpPion")->GetValue();
                    trackFromTPCPids->dedxexpkaon = (*tr)->Get< COMET::IRealDatum >("tpcPid_dEdxexpKaon")->GetValue();
                    trackFromTPCPids->dedxexpproton = (*tr)->Get< COMET::IRealDatum >("tpcPid_dEdxexpProton")->GetValue();
                    trackFromTPCPids->sigmaele = (*tr)->Get< COMET::IRealDatum >("tpcPid_SigmaEle")->GetValue();
                    trackFromTPCPids->sigmamuon = (*tr)->Get< COMET::IRealDatum >("tpcPid_SigmaMuon")->GetValue();
                    trackFromTPCPids->sigmapion = (*tr)->Get< COMET::IRealDatum >("tpcPid_SigmaPion")->GetValue();
                    trackFromTPCPids->sigmakaon = (*tr)->Get< COMET::IRealDatum >("tpcPid_SigmaKaon")->GetValue();
                    trackFromTPCPids->sigmaproton = (*tr)->Get< COMET::IRealDatum >("tpcPid_SigmaProton")->GetValue();
                    trackFromTPCPids->pullele = (*tr)->Get< COMET::IRealDatum >("tpcPid_PullEle")->GetValue();
                    trackFromTPCPids->pullmuon = (*tr)->Get< COMET::IRealDatum >("tpcPid_PullMuon")->GetValue();
                    trackFromTPCPids->pullpion = (*tr)->Get< COMET::IRealDatum >("tpcPid_PullPion")->GetValue();
                    trackFromTPCPids->pullkaon = (*tr)->Get< COMET::IRealDatum >("tpcPid_PullKaon")->GetValue();
                    trackFromTPCPids->pullproton = (*tr)->Get< COMET::IRealDatum >("tpcPid_PullProton")->GetValue();

                    //dEdx->Fill((*tr)->Get< COMET::IRealDatum >("tpcPid_Ccorr")->GetValue());
                }

                COMET::IHandle<COMET::ITrackState> frontstate ((*(nodes.begin()))->GetState());

                COMET::IHandle<COMET::ITrackState> backstate ((*(nodes.rbegin()))->GetState());

                double txin,tyin,etxin,etyin;
                double txfin,tyfin,etxfin,etyfin;
                getTangents(frontstate,txin,etxin,tyin,etyin);
                getTangents(backstate,txfin,etxfin,tyfin,etyfin);

                trackFromTPCPids->txin=txin;
                trackFromTPCPids->tyin=tyin;
                trackFromTPCPids->txfin=txfin;
                trackFromTPCPids->tyfin=tyfin;

                trackFromTPCPids->npoint=nodes.size();

                int l=0;

                TVector3 dir = track->GetDirection();
                TLorentzVector pos = track->GetPosition();

                double T0 = pos.T() + fTimeOffset;

                //COMETDebug(T0);

                trackFromTPCPids->T0 = T0;
                trackFromTPCPids->theta = dir.Theta();
                trackFromTPCPids->phi = dir.Phi();

                TVector3 dirfront = frontstate->GetDirection();
                TLorentzVector posfront = frontstate->GetPosition();

                double xf,yf,zf,rhofront,tanphi,txfront,tyfront;
                xf = posfront[0];
                yf = posfront[1];
                zf = posfront[2];
                rhofront =  1./frontstate->GetCurvature();
                tanphi =  dirfront[1]/dirfront[2];
                txfront = dirfront[0]/dirfront[2];
                tyfront = dirfront[1]/dirfront[2];
                double phi = TMath::ATan(tanphi);

                //#ifdef MC
                if(event.GetContext().IsMC()) {
                    //if( nodes.size() < 35 ) continue;
                    COMET::IG4Trajectory g4traj;

                    COMET::IG4TrajectoryPoint G4Point = GetG4Point(event,xf,yf,zf,rhofront,phi,txfront,g4traj);

                    COMET::IG4Trajectory G4trackTPC=GetTrajectory(event, *track);

                    TLorentzVector trinitpos = G4trackTPC.GetInitialPosition();
                    TLorentzVector trfinpos = G4trackTPC.GetInitialPosition();
                    //Int_t truecharge = G4trackTPC.GetParticle()->Charge();

                    //trajpvect = trajpair.second.GetTrajectoryPoints();

                    // Unused variable therefore commenting out to remove compiler warnings
                    //int parentid = g4traj.GetParentId();
                    //int parid = g4traj.GetTrackId();
                    int parcode = g4traj.GetPDGEncoding();

                    //GetTrueVertex(event);

                    //COMETDebug("track initial position"<<trinitpos.X()<<" "<<trinitpos.Y()<<" "<<trinitpos.Z());
                    //COMETDebug("track info "<<parentid<<" "<<parid<<" "<<parcode);
                    //COMETDebug("MC PID "<<G4trackTPC.GetPDGEncoding());
                    //char s[1];
                    //gets(s);

                    trackFromTPCPids->TruePID=parcode;
                    trackFromTPCPids->TrueX=trinitpos.X();
                    trackFromTPCPids->TrueY=trinitpos.Y();
                    trackFromTPCPids->TrueZ=trinitpos.Z();
                    trackFromTPCPids->TrueP=G4Point.GetMomentum().Mag();
                }
                //#endif

                if( dir.Theta() > 3.141592/2.-0.01 ) {

                    //COMETDebug( " THETA " );
                    //COMETDebug( " N Event " << nevt );
                    //COMETDebug( " X " << pos.X() );
                    //COMETDebug( " Y " << pos.Y() );
                    //COMETDebug( " Z " << pos.Z() );
                    //COMETDebug( " Dx" <<  dir.X() );
                    //COMETDebug( " Dy" <<  dir.Y() );
                    //COMETDebug( " Dz" <<  dir.Z() );
                }

                //theta->Fill(dir.Theta());
                //thetavsz->Fill(pos.Z(),dir.Theta());
                //phiangle->Fill(dir.Phi());
                //X->Fill(pos.X());
                //Y->Fill(pos.Y());
                //Z->Fill(pos.Z());
                //NHits->Fill((double)nodes.size());
                //Goodness->Fill(track->GetQuality());

                bool flag = false;

                for( COMET::IReconNodeContainer::iterator  ni = nodes.begin(); ni < nodes.end();ni++) {

                    flag = false;

                    COMET::IHandle<COMET::ITrackState> state = (*ni)->GetState();
                    COMET::IHandle<COMET::IReconBase>  rbase = (*ni)->GetObject();
                    //double charge=state->GetEDeposit();

                    getTangents(state, tx, etx, ty, ety);

                    COMET::IHandle<COMET::IHitSelection> hits = rbase->GetHits();

                    COMET::ITPCGeom& tpcGeom = COMET::IGeomInfo::TPC();

                    COMET::IHitSelection::const_iterator cHit = (*hits).begin();

                    COMET::IHandle<COMET::IComboHit> combo(*cHit);

                    int count = 0;

                    for(COMET::IHitSelection::const_iterator Hit = combo->GetHits().begin() ;Hit!=combo->GetHits().end();Hit++) {
                                    ///check if the cluster is on the board of 1 MM module or on the cathod
                        count++;

                        //int node=(*Hit)->GetGeomId();
                        //int key=(*Hit)->GetGeoKey();

                        int padId = tpcGeom.GeomIdToPad((*Hit)->GetGeomId());
                        //int padId = tpcGeom.ChannelToPad(node,key);

                        //int MMmod = tpcGeom.GeomIdToMM((*Hit)->GetGeomId());
                        //int TPCHalf = tpcGeom.GeomIdToHalf((*Hit)->GetGeomId());
                        //int TPC = tpcGeom.GeomIdToTPC((*Hit)->GetGeomId());

                        TVector3 padpos;

                        tpcGeom.GeomIdToGlobalXYZ((*Hit)->GetGeomId(),padpos);

                        //int modId = (gGeoManager->FindNode(padpos.X(),padpos.Y(),padpos.Z()))->GetNumber();

                        //int column = tpcGeom.PadToColumn(padId);
                        int row = tpcGeom.PadToRow(padId);

                        if(row==0 || row==47)
                            flag=true;

                        //COMETDebug("count "<<count<<" X "<<padpos.X()<<" "<<MMmod<<" "<<TPCHalf<<" "<<TPC<<" column "<<column<<" row "<<row);

                    }

                    //Charge->Fill( state->GetEDeposit() );

                    //if( ni == nodes.begin() ) Curvature->Fill(track->GetCurvature());

                    l++;

                }

                if( !frontstate || !backstate ) continue;

                double momentum=0;
                double errmomentum=0;

                getMomentum(frontstate,momentum,errmomentum);
                //p = (0.3*unit::GeV/unit::tesla/unit::m)*Bmag/(frontstate->GetCurvature()*TMath::Sqrt(1.-frontstate->GetDirection()[0]*frontstate->GetDirection()[0]));
                trackFromTPCPids->curvature = track->GetCurvature();
                trackFromTPCPids->momentum = momentum;
                trackFromTPCPids->momentumerr = errmomentum;

                TVector3 frontdir = frontstate->GetDirection();
                TLorentzVector frontpos = frontstate->GetPosition();
                TVector3 backdir = backstate->GetDirection();
                TLorentzVector backpos = backstate->GetPosition();

                //q = ( frontstate->GetCurvature()>0? -1:1 );
                x = frontpos[0];
                y = frontpos[1];
                z = frontpos[2];
                //ty = frontdir[1]/frontdir[2];
                //tx = frontdir[0]/frontdir[2];
                //rho = frontstate->GetCurvature();
                //p = 1./fTPCRecon->GetTPCLikFit()->Rho2invp(TMath::Abs(rho),tx,ty);

                trackFromTPCPids->zin=z;
                trackFromTPCPids->xin=x;
                trackFromTPCPids->yin=y;

                trackFromTPCPids->zfin=backpos[2];
                trackFromTPCPids->yfin=backpos[1];
                trackFromTPCPids->xfin=backpos[0];

                //COMETDebug("Event ID "<<evtID<<" CT meas "<<trackFromTPCPids->dedxcorr[itr]<<" momentum "<<tree.momentum[itr]<<" err momentum "<<tree.momentumerr[itr]<<" x in "<<tree.xin[itr]<<" x fin "<<tree.xfin);

                itr++;
                fNTracksFromTPCPids++;

            }                       //tracks loop

        }                           //end loop if track exist
    }                               //end loop if FitResult
    else {
        //COMETDebug( " No Fit " );
    }

    //treefile->Fill();

    return true;
}


COMET::IExample2010aAnalysis1Module::ITrackFromTPCPids::ITrackFromTPCPids(){ }
COMET::IExample2010aAnalysis1Module::ITrackFromTPCPids::~ITrackFromTPCPids(){ }

COMET::IExample2010aAnalysis1Module::IFGDHit::IFGDHit(){ }
COMET::IExample2010aAnalysis1Module::IFGDHit::~IFGDHit(){ }

COMET::IExample2010aAnalysis1Module::IP0DHit::IP0DHit(){ }
COMET::IExample2010aAnalysis1Module::IP0DHit::~IP0DHit(){ }

COMET::IExample2010aAnalysis1Module::ISMRDHit::ISMRDHit(){ }
COMET::IExample2010aAnalysis1Module::ISMRDHit::~ISMRDHit(){ }

COMET::IExample2010aAnalysis1Module::ITrueNuVertex::ITrueNuVertex(){ }
COMET::IExample2010aAnalysis1Module::ITrueNuVertex::~ITrueNuVertex(){ }

void COMET::IExample2010aAnalysis1Module::getTangents( COMET::IHandle<COMET::ITrackState> state,double &tx,double &etx,double &ty,double &ety) {

    tx =  state->GetDirection()[0]/state->GetDirection()[2];
    ty =  state->GetDirection()[1]/state->GetDirection()[2];

    CLHEP::HepMatrix  Jacobian(2,3,0);
    CLHEP::HepMatrix  Covariance(3,3,0);
    CLHEP::HepMatrix  Cov(2,2,0);

    for(int i = 0; i < 3; i++ )
    for(int j = 0; j < 3; j++ ) {
        Covariance[i][j] = state->GetCovarianceValue(state->GetDirectionIndex()+i,state->GetDirectionIndex()+j);
    }

    Jacobian[0][0] = 1./state->GetDirection()[2];
    Jacobian[0][1] = 0.;
    Jacobian[0][2] = -tx/state->GetDirection()[2];

    Jacobian[1][0] = 0.;
    Jacobian[1][1] = 1./state->GetDirection()[2];
    Jacobian[1][2] = -ty/state->GetDirection()[2];

    Cov = Jacobian*Covariance*Jacobian.T();

    etx = TMath::Sqrt(Cov[0][0]);
    ety = TMath::Sqrt(Cov[1][1]);

    return;

}


void COMET::IExample2010aAnalysis1Module::getMomentum( COMET::IHandle<COMET::ITrackState> state,double &p,double &ep) {

    COMET::IHandle<COMET::IReconState> bstate = state;
    p  = TrackingUtils::GetMomentum(bstate);
    ep = TrackingUtils::GetMomentumError(bstate);

    return;
}


COMET::IG4Trajectory COMET::IExample2010aAnalysis1Module::GetTrajectory(COMET::ICOMETEvent& event,COMET::IReconTrack &track) {
    //COMETDebug("****** GET ASSOCIATE TRAJECTORY TO THE TRACK ****** ");
    COMET::IG4Trajectory traj;
    COMET::IHandle<COMET::IHitSelection> hits = track.GetHits();
    COMET::IHandle<COMET::IG4TrajectoryContainer> trajectories = event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");
    if (!trajectories) return COMET::IG4Trajectory();
    std::vector<int> search;
    search.resize(trajectories->size());
    if(hits->size()>0) {
        // Loop over reconstructed hits
        for( COMET::IHitSelection::const_iterator ht = hits->begin(); ht != hits->end(); ht++ ) {
            //            const COMET::IHandle<COMET::IHit> hit=ht;
            COMET::IHandle<COMET::IComboHit> ch(*ht);
            if(!ch)
                continue;
            // Loop over hits in Combo
            for (COMET::IHitSelection::const_iterator ch_member = ch->GetHits().begin();
            ch_member != ch->GetHits().end(); ++ch_member) {

                COMET::IMultiHit *mh = dynamic_cast<COMET::IMultiHit*> (&(**ch_member));

                if( !mh ) {         // this is the old method with no waveforms.

                    // Get monte Carlo hit associated.
                    COMET::IHandle<COMET::IMCHit> mch(*ch_member);
                    if (!mch) {
                        continue;
                    }
                    for (std::vector<COMET::IG4VHit*>::const_iterator h = (mch->GetContributors()).begin(); h != (mch->GetContributors()).end(); ++h) {
                        if( !*h ) continue;
                        COMET::IG4HitSegment *g4hitseg = dynamic_cast<COMET::IG4HitSegment*> (*h);
                        if( g4hitseg != NULL ) {
                            for( unsigned int i = 0; i != g4hitseg->GetContributors().size(); ++i) {
                                unsigned int vv = g4hitseg->GetContributors()[i];
                                if( vv >= search.size() )
                                    search.resize(vv+1);
                                search[vv]=search[vv]+1;
                            }
                        }
                        else
                            COMETWarn(  " NO HITS " );
                    }
                }                   ///end old method

                else {              // This is the new method with waveform
                    typedef std::vector< COMET::IHandle< COMET::ISingleHit > > Hits;
                    typedef Hits::const_iterator iterator;

                    for( iterator wfit = mh->begin() ; wfit != mh->end(); wfit++ ) {

                                    // Get monte Carlo hit associated.
                        COMET::IHandle<COMET::IMCHit> mch(*wfit);

                        if( mch == NULL ) continue;

                        for (std::vector<COMET::IG4VHit*>::const_iterator h = (mch->GetContributors()).begin();
                        h != (mch->GetContributors()).end(); ++h) {

                            if( !*h ) continue;

                            COMET::IG4HitSegment *g4hitseg = dynamic_cast<COMET::IG4HitSegment*> (*h);

                            if( g4hitseg != NULL ) {

                                for( unsigned int i = 0; i != g4hitseg->GetContributors().size(); ++i) {
                                    unsigned int vv = g4hitseg->GetContributors()[i];

                                    if( vv >= search.size() )
                                        search.resize(vv+1);

                                    search[vv]=search[vv]+1;
                                }
                            }
                            else
                                COMETWarn(  " NO HITS " );
                        }
                    }

                }                   //end of the new method

            }                       //close some loops
        }
    }

    int maxim = 0;
    int imax = -1;
    for( unsigned int i = 0; i < search.size() ; i++ ) {

        //COMETDebug("searcha lla fine "<<i<<" search "<<search[i]<<" maxim "<<maxim);

        if( search[i] > maxim ) {
            maxim = search[i];
            imax = i;
        }
    }

    //COMETDebug("imax "<<imax);

    if( imax >= 0 ) {
        COMETDebug(  "...... trajectory found " );
        return (*trajectories)[imax];
    }
    else {
        COMETWarn("...... trajectory not found for this track " << hits->size() << "  " <<  search.size());
        return COMET::IG4Trajectory();
    }

    return COMET::IG4Trajectory();
}


COMET::IG4TrajectoryPoint COMET::IExample2010aAnalysis1Module::GetG4Point(COMET::ICOMETEvent& event,double x0,double y0,double z0,double rho, double phi, double tx,COMET::IG4Trajectory &g4traj) {

    double mindist = 100000000.;
    COMET::IG4TrajectoryContainer fTrajCon  =  *(event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories"));
    Int_t parid, parentid, parcode;
    std::pair<Int_t,COMET::IG4Trajectory> trajpair;
    std::vector<COMET::IG4TrajectoryPoint> trajpvect;
    COMET::IG4TrajectoryPoint g4Point;

    for (COMET::IG4TrajectoryContainer::const_iterator t = fTrajCon.begin(); t != fTrajCon.end(); t++) {
        trajpair = *t;
        trajpvect = trajpair.second.GetTrajectoryPoints();
        parentid = trajpair.second.GetParentId();
        parid = trajpair.second.GetTrackId();
        parcode = trajpair.second.GetPDGEncoding();

        if( trajpair.second.GetParticle() == NULL ) continue;

        if(  trajpair.second.GetParticle()->Charge() == 0 ) continue;

        for( std::vector<COMET::IG4TrajectoryPoint>::iterator tp = trajpvect.begin(); tp != trajpvect.end(); tp++ ) {

            //     double x = x0 + tx*(tp->GetPosition().Z()-z0);
            double y = y0 + 1./rho * ( TMath::Cos( TMath::ASin(-(tp->GetPosition().Z()-z0)*rho + TMath::Sin(phi)) )-TMath::Cos(phi) );
            // double z = z0;

            double dist = TMath::Sqrt( pow( TMath::Abs(tp->GetPosition().Y()-y), 2.) + pow(TMath::Abs(tp->GetPosition().Z()-z0), 2.));

            if(  dist  < mindist ) {
                mindist = dist;
                g4Point = *tp;
                g4traj = trajpair.second;

            }
        }
    }

    return g4Point;
}


int COMET::IExample2010aAnalysis1Module::GetNeutrinoFamily (std::string fReactionName) {

    int nuFamily = -1;
    size_t foundProcess;

    //Charged current
    //------------------------------------------------
    //Check if CC_QEL
    foundProcess = fReactionName.find("nu:14");
    if(foundProcess != std::string::npos)
        nuFamily = 14;

    foundProcess = fReactionName.find("nu:-14");
    if(foundProcess != std::string::npos)
        nuFamily = -14;

    foundProcess = fReactionName.find("nu:12");
    if(foundProcess != std::string::npos)
        nuFamily = 12;

    foundProcess = fReactionName.find("nu:-12");
    if(foundProcess != std::string::npos)
        nuFamily = -12;

    return nuFamily;
}


int COMET::IExample2010aAnalysis1Module::GetReactionCode (std::string fReactionName) {

    int reactionCode = -1;
    size_t foundProcess;

    //Charged current
    //------------------------------------------------
    //Check if CC_QEL
    foundProcess = fReactionName.find("Weak[CC],QES");
    if ( foundProcess != std::string::npos )
        reactionCode = 1;

    //Check if CC_RES
    foundProcess = fReactionName.find("Weak[CC],RES");
    if ( foundProcess != std::string::npos )
        reactionCode = 2;

    //Check if CC_DIS
    foundProcess = fReactionName.find("Weak[CC],DIS");
    if ( foundProcess != std::string::npos )
        reactionCode = 3;

    //Check if CC_COHpi
    foundProcess = fReactionName.find("Weak[CC],COHPi");
    if ( foundProcess != std::string::npos )
        reactionCode = 4;

    //Check if CC_COHel
    foundProcess = fReactionName.find("Weak[CC],COHEl");
    if ( foundProcess != std::string::npos )
        reactionCode = 5;

    //Neutral current
    //------------------------------------------------
    //Check if NC_QEL
    foundProcess = fReactionName.find("Weak[NC],QES");
    if ( foundProcess != std::string::npos )
        reactionCode = 11;

    //Check if NC_RES
    foundProcess = fReactionName.find("Weak[NC],RES");
    if ( foundProcess != std::string::npos )
        reactionCode = 12;

    //Check if NC_DIS
    foundProcess = fReactionName.find("Weak[NC],DIS");
    if ( foundProcess != std::string::npos )
        reactionCode = 13;

    //Check if NC_COHpi
    foundProcess = fReactionName.find("Weak[NC],COHPi");
    if ( foundProcess != std::string::npos )
        reactionCode = 14;

    //Check if NC_COHel
    foundProcess = fReactionName.find("Weak[NC],COHEl");
    if ( foundProcess != std::string::npos )
        reactionCode = 15;

    //Other
    //------------------------------------------------
    //Check if AMNuGamma
    foundProcess = fReactionName.find("AMNuGamma");
    if ( foundProcess != std::string::npos )
        reactionCode = 21;

    //Check if NuEEL
    foundProcess = fReactionName.find("NuEEL");
    if ( foundProcess != std::string::npos )
        reactionCode = 22;

    //Check if IMD
    foundProcess = fReactionName.find("IMD");
    if ( foundProcess != std::string::npos )
        reactionCode = 23;

    //Charm production
    //------------------------------------------------
    foundProcess = fReactionName.find("charm");
    if ( foundProcess != std::string::npos )
        reactionCode = 31;

    foundProcess = fReactionName.find("QES;charm");
    if ( foundProcess != std::string::npos )
        reactionCode = 32;

    foundProcess = fReactionName.find("DIS;charm");
    if ( foundProcess != std::string::npos )
        reactionCode = 33;

    return reactionCode;
}


void COMET::IExample2010aAnalysis1Module::GetTrueVertex(COMET::ICOMETEvent& event) {

    Int_t nu=0;

    COMET::IHandle<COMET::IG4PrimaryVertexContainer> g4PrimVertex00 = event.Get<COMET::IG4PrimaryVertexContainer>("truth/G4PrimVertex00");

    COMETDebug("primary vertex "<<g4PrimVertex00->size());

    for(std::vector<COMET::IG4PrimaryVertex>::const_iterator vertexIter = g4PrimVertex00->begin(); vertexIter != g4PrimVertex00->end(); ++vertexIter) {

        TLorentzVector vertexPosition = vertexIter->GetPosition();
        //vtxToFill->Vertex = vertexIter->GetPosition();
        std::string Generator = vertexIter->GetGeneratorName();

        std::string ReactionCode = vertexIter->GetReaction();

        int codeneut=std::atoi(ReactionCode.c_str());

        //TLorentzVector nuene = vertexIter->GetEnergy();
        //double NeutrinoEnergy = sqrt(pow(nuene.X(),2)+pow(nuene.Y(),2)+pow(nuene.Z(),2));

        //COMETDebug(ReactionCode);
        //COMETDebug(" Z "<<vertexIter->GetPosition().Z()<<" Y "<<vertexIter->GetPosition().Y()<<" X "<<vertexIter->GetPosition().X());

        COMET::TG4PrimaryParticleContainer primaryParticlesCopy(vertexIter->GetPrimaryParticles());

        //sort(primaryParticlesCopy.begin(),primaryParticlesCopy.end(),DescendingKE);

        for(std::vector<COMET::IG4PrimaryVertex>::const_iterator infoVtxIter = vertexIter->GetInfoVertex().begin();
        infoVtxIter != vertexIter->GetInfoVertex().end(); ++infoVtxIter) {

            for(std::vector<COMET::IG4PrimaryParticle>::const_iterator particleIter = infoVtxIter->GetPrimaryParticles().begin();
                particleIter != infoVtxIter->GetPrimaryParticles().end();
            ++particleIter) {

                // particleIter is a pointer to a COMET::IG4PrimaryParticle
                Int_t pdg = particleIter->GetPDGCode();

                if(fabs(pdg)==12 || fabs(pdg)==14 || fabs(pdg)==16) {
                    TLorentzVector nuene = particleIter->GetMomentum();
                    Float_t NeutrinoEnergy = sqrt(pow(nuene.X(),2)+pow(nuene.Y(),2)+pow(nuene.Z(),2));

                    int FlorCode, NuCode;

                    /* if(GENIE)
                      {

                        FlorCode = GetReactionCode(ReactionCode);
                        NuCode = GetNeutrinoFamily(ReactionCode);

                        COMETDebug("reaction type "<<FlorCode<<" nu family "<<NuCode);

                      }
                    else
                    */
                    {
                        FlorCode = codeneut;
                        NuCode = pdg;
                    }

                    COMETDebug("found a neutrino "<<nu);
                    COMETDebug(ReactionCode<<" "<<FlorCode<<" "<<NuCode);
                    COMETDebug("pdg "<<pdg<<" Z "<<vertexIter->GetPosition().Z());
                    COMETDebug("True neutrino energy "<<particleIter->GetMomentum().Z()<<" "<<NeutrinoEnergy);

                    // Create neutrino vertex object to fill
                    ITrueNuVertex* trueNuVertex = new((*fTrueNuVertices)[fNTrueNuVertices]) ITrueNuVertex;

                    trueNuVertex->TrueVertexX=vertexIter->GetPosition().X();
                    trueNuVertex->TrueVertexY=vertexIter->GetPosition().Y();
                    trueNuVertex->TrueVertexZ=vertexIter->GetPosition().Z();
                    trueNuVertex->TrueNuEnergy=NeutrinoEnergy;
                    trueNuVertex->TrueReaction=FlorCode;
                    trueNuVertex->TrueFamily=NuCode;
                    fNTrueNuVertices++;

                }

            }
        }
    }

}


void COMET::IExample2010aAnalysis1Module::FillFGDHits(COMET::IHandle<COMET::IHitSelection> fgd) {
    //**************************************************************

    //------------ Draw reconstructed FGD hits ------------------------------
    if( !fgd ) return;

    for(COMET::IHitSelection::const_iterator hsel = fgd->begin(); hsel != fgd->end(); hsel++) {
        // Create FGD hit object to fill
        IFGDHit* fgdHit = new((*fFGDHits)[fNFGDHits]) IFGDHit;

        const TVector3 pos = (*hsel)->GetPosition();

        fgdHit->fgdCharge=(*hsel)->GetCharge();

        fgdHit->fgdT=(*hsel)->GetTime();

        //COMETDebug(fgdHit->fgdT);

        if( (*hsel)->IsYHit() ) {

            fgdHit->fgdX=0;
            fgdHit->fgdY=pos.Y();
            fgdHit->fgdZ=pos.Z();
        }
        else if( (*hsel)->IsXHit()     ) {
            fgdHit->fgdX=pos.X();
            fgdHit->fgdY=0;
            fgdHit->fgdZ=pos.Z();

        }

        //COMETDebug(countfgd<<" "<<fgdHit->fgdZ[countfgd]<<" "<<tree.fgdX[countfgd]<<" "<<tree.fgdY[countfgd]<<" "<<tree.fgdCharge);

        fNFGDHits++;
        //char s[1];
        //gets(s);

    }

}


void COMET::IExample2010aAnalysis1Module::FillP0DHits(COMET::IHandle<COMET::IHitSelection> p0d) {
    //**************************************************************

    //------------ Draw reconstructed P0D hits ------------------------------
    if( !p0d ) {
        return;
    }

    // tree.p0dSize=p0d->size();

    for(COMET::IHitSelection::const_iterator hsel = p0d->begin(); hsel != p0d->end(); hsel++) {

        // Create P0D hit object to fill
        IP0DHit* p0dHit = new((*fP0DHits)[fNP0DHits]) IP0DHit;

        if((*hsel)->GetTime()<30000) {

            const TVector3 pos = (*hsel)->GetPosition();

            p0dHit->p0dCharge=(*hsel)->GetCharge();

            p0dHit->p0dT=(*hsel)->GetTime();

            //COMETDebug(p0dHit->p0dT);

            if( (*hsel)->IsYHit() ) {

                p0dHit->p0dX=0;
                p0dHit->p0dY=pos.Y();
                p0dHit->p0dZ=pos.Z();
            }
            else if( (*hsel)->IsXHit()     ) {
                p0dHit->p0dX=pos.X();
                p0dHit->p0dY=0;
                p0dHit->p0dZ=pos.Z();

            }

            //COMETDebug(countp0d<<" "<<p0dHit->p0dZ[countp0d]<<" "<<tree.p0dX[countp0d]<<" "<<tree.p0dY[countp0d]<<" "<<tree.p0dCharge);

            fNP0DHits++;
            //char s[1];
            //gets(s);
        }

    }

}


//**************************************************************

void COMET::IExample2010aAnalysis1Module::FillSMRDHits(COMET::IHandle<COMET::IHitSelection> smrd) {
    //**************************************************************
    //------------ Fill reconstructed SMRDD hits ----------------
    if( !smrd ) {
        return;
    }

    for(COMET::IHitSelection::const_iterator hsel = smrd->begin(); hsel != smrd->end(); hsel++) {

        // Create SMRD hit object to fill
        ISMRDHit* smrdHit = new((*fSMRDHits)[fNSMRDHits]) ISMRDHit;

        if((*hsel)->GetTime()<30000) {

            const TVector3 pos = (*hsel)->GetPosition();

            smrdHit->smrdCharge=(*hsel)->GetCharge();

            smrdHit->smrdT=(*hsel)->GetTime();

            /*     COMETDebug(smrdHit->smrdCharge); */

            if( (*hsel)->IsYHit() ) {

                smrdHit->smrdX=0;
                smrdHit->smrdY=pos.Y();
                smrdHit->smrdZ=pos.Z();
            }
            else if( (*hsel)->IsXHit()     ) {
                smrdHit->smrdX=pos.X();
                smrdHit->smrdY=0;
                smrdHit->smrdZ=pos.Z();

            }

            //COMETDebug(countsmrd<<" "<<smrdHit->smrdZ[countsmrd]<<" "<<tree.smrdX[countsmrd]<<" "<<tree.smrdY[countsmrd]<<" "<<tree.smrdCharge);

            fNSMRDHits++;
            //char s[1];
            //gets(s);
        }

    }

}
#endif // OAANALYSIS_WITHOUT_RECON_CALIB
// set the environment variable 
// CMTEXTRATAGS to OAANALYSIS_WITHOUT_RECON_CALIB
// to avoid this code