GetClusterFitExtrapolatedTOFPixelID.hpp

//__________________________________________________________________________________________________________________________
std::vector<int>
GetClusterFitExtrapolatedTOFPixelID(const Track& track,
                                    const RecEvent& recEvent,
                                    const det::TOFWall& TOFwall,
                                    bool debug = false)
{
    if (debug == true)
        cerr<<"\n\nNumber of points on track in the MTPC = "<<track.GetNumberOfClusters(TrackConst::eMTPC)<<endl;

    const utl::Point& last_point_on_track = track.GetLastPointOnTrack();
    std::vector<Point> clusters;


    int iPoint = 0;
    for (ClusterIndexIterator ClusterIter = track.ClustersEnd(); ClusterIter != track.ClustersBegin(); --ClusterIter)
    {
        if (clusters.size() > 9)
            break;

        if (debug == true)
            cerr<<"\nClusterIter =  "<<*ClusterIter;

        if (iPoint == 0) --ClusterIter;
            ++iPoint;

        if (debug == true)
            cerr<<"\nFor iPoint =  "<<iPoint;

        const Cluster& cluster = recEvent.Get(*ClusterIter);
        clusters.push_back(cluster.GetPosition());

        if (debug == true)
        {
            cerr<<":  "<<cluster.GetPosition().GetX()<<", "<<cluster.GetPosition().GetY()<<", "<<cluster.GetPosition().GetZ();

            if (last_point_on_track == cluster.GetPosition())
                cerr<<".  Found last point on track, corresponding to iPoint = "<<iPoint;
        }
    }

    if (debug == true) cerr<<"\n";


    /*
    int pointsIter = 0;
    for (ClusterIndexIterator ClusterIter = track.ClustersEnd(); pointsIter < 10; ++pointsIter)
    {
        if (debug == true)
            cerr<<"\nClusterIter =  "<<*ClusterIter;

        const Cluster& cluster = recEvent.Get(*(--ClusterIter));
        clusters.push_back(cluster.GetPosition());

        if (debug == true)
        {
            cerr<<"\nFor pointsIter =  "<<pointsIter;
            cerr<<":  "<<cluster.GetPosition().GetX()<<", "<<cluster.GetPosition().GetY()<<", "<<cluster.GetPosition().GetZ();

            if (last_point_on_track == cluster.GetPosition())
                cerr<<".  Found last point on track, corresponding to pointsIter = "<<pointsIter;
        }
    }

    if (debug == true) cerr<<"\n";
    */


    TGraph * graph_track_XZ = new TGraph(clusters.size());
    TGraph * graph_track_YZ = new TGraph(clusters.size());
    for(unsigned int iclust = 0; iclust < clusters.size(); iclust++)
    {
        graph_track_XZ->SetPoint((int)iclust, clusters.at(iclust).GetZ(), clusters.at(iclust).GetX());
        graph_track_YZ->SetPoint((int)iclust, clusters.at(iclust).GetZ(), clusters.at(iclust).GetY());
    }
    TF1 * lin_fit_XZ = new TF1("track_XZ", "[0]*x+[1]", 650, 830);
    TF1 * lin_fit_YZ = new TF1("track_XZ", "[0]*x+[1]", 650, 830);


    graph_track_XZ->Fit(lin_fit_XZ, "Q"); //"V"
    graph_track_YZ->Fit(lin_fit_YZ, "Q");

    // TFitResultPtr r1 = graph_track_XZ->Fit(lin_fit_XZ, "SQ");
    // TFitResultPtr r2 = graph_track_YZ->Fit(lin_fit_YZ, "SQ");
    // cout<<"Chi-sq/ndf:  Fit XZ = "<<r1->Chi2()/r1->Ndf()<<", Fit YZ = "<<r2->Chi2()/r2->Ndf()<<endl;
    // delete r1;
    // delete r2;


    std::vector<int> vTOFHitIndex;
    for(int hitindex = 1; hitindex <= 891; ++hitindex)
    {
        const det::TOFScintillator& Iscint = TOFwall.GetScintillator(hitindex);
        double extrapolationPositionZ = Iscint.GetCenterPosition().GetZ() - 2.3/2;
        double extrapolationPositionX = lin_fit_XZ->Eval(extrapolationPositionZ);
        double extrapolationPositionY = lin_fit_YZ->Eval(extrapolationPositionZ);

        if (extrapolationPositionX < Iscint.GetCenterPosition().GetX() - Iscint.GetWidth()/2.) continue;
        if (extrapolationPositionX > Iscint.GetCenterPosition().GetX() + Iscint.GetWidth()/2.) continue;
        if (extrapolationPositionY < Iscint.GetCenterPosition().GetY() - 3.4/2) continue;
        if (extrapolationPositionY > Iscint.GetCenterPosition().GetY() + 3.4/2) continue;

        vTOFHitIndex.push_back(hitindex);
    }


    delete graph_track_XZ;
    delete graph_track_YZ;
    delete lin_fit_XZ;
    delete lin_fit_YZ;

    return vTOFHitIndex;
}