tests/src/ExtractEfficiencies.cpp

Go to the documentation of this file.

#include "MapsTrackManager.hh"
#include "MapsTrack.hh"
#include "MapsSensor.hh"
#include "MapsException.hh"
#include <cmath>
#include <algorithm>
#include <functional>
#include <iostream>
#include "TRandom2.h"
#include "TFile.h"
#include "TH1F.h"
#include "Operators.h"

/** ExtractEfficiencies.cpp
 * 
 * Workhorse for determining sensor efficiency once tracks have been processed
 * and you're happy with the alignment.
 * 
 * Exports plots to supplied outpufile name.
 * 
 * This is a good place to start to understand the general code structure, and
 * as a general example for how you could write your own tool.
 * 
 * Usage: ExtractEfficiencies input file outputfile
 * 
 * Jamie Ballin, Imperial College London
 *              February 2008
 * 
 */

int main(int argc, const char **argv) {

        if (argc != 3) {
                std::cout << "Usage: ExtractEfficiencies <input file> <outputfile>\n";
                return 0;
        }
        std::cout << "Welcome to ExtractEfficiencies...\n";

        //Define chi squared probability cuts
        double cut(0.05);
        //Define errors for chi squared calculation
        std::pair<double, double> error(0.018, 0.018);
        //std::pair<double, double> error(0.026, 0.02);
        std::cout << "\tUsing cuts for probability of: \t" << cut << "\n";
        std::cout << "\tError parameters: \t" << error << "\n";

        std::cout << "\tRecreating from root file...\n";
        char input[100];
        strcpy(input, argv[1]);

        //Create a maps track manager to recreate the tracks
        MapsTrackManager mtm2;
        //Reinitialise all the tracks and sensors in mtm2
        mtm2.recreateFromRootFile(input);

        std::vector<MapsTrack*>& tracks = mtm2.getTracks();
        std::vector<MapsSensor*>& sensors = mtm2.getSensors();
        MapsSensor* s1;
        MapsSensor* s2;
        for (std::vector<MapsSensor*>::const_iterator it = sensors.begin(); it
                        != sensors.end(); ++it) {
                if ((*it)->id() == 8) {
                        s1 = *it;
                } else if ((*it)->id() == 6) {
                        s2 = *it;
                }
        }

        //count efficient and inefficient frequency
        unsigned ineff(0), eff(0), deadArs(0);
        //counts specific to shapers and samplers
        unsigned shapEff(0), sampEff(0);
        unsigned shapIneff(0), sampIneff(0);
        unsigned candidateTracks(0);

        for (std::vector<MapsTrack*>::iterator it = tracks.begin(); it
                        != tracks.end(); it++) {
                MapsTrack* t = *it;

                //this is not obligatory, but highly advised!
                t->setTrackError(error);

                //see if we have a good three hit track first: for tracks which are just three hits anyway
                //t3 will be a copy of t
                MapsTrack t3;
                //construct the three hit track
                t->make3HitTrack(t3);
                t3.setTrackError(error);

                //evaluate quality of the 3 hit track
                if (t3.chiSqProb(0) > cut && t3.chiSqProb(1) > cut) {
                        //it's a real track then
                        candidateTracks++;

                        //get the original track's fourth sensor
                        MapsSensor* s4 = t->fourthSensor();

                        //pose the ultimate question...
                        std::pair<double, double> resid;
                        unsigned fourthConfirm = t->getFourthHitResidual(t3, resid);

                        if (fourthConfirm == 0 && t->chiSqProb(0) > cut && t->chiSqProb(1)
                                        > cut) {
                                //fourth sensor was efficient
                                //one could also cut on the residual value, if you wanted
                                s4->registerTrackConfirm(t->fourthSensorThresh(), t->getHits().at(s4), t->timeStamp(), resid);
                                eff++;
                                if (t->getHits().at(s4).first < 84)
                                        shapEff++;
                                else
                                        sampEff++;
                        } else {
                                //fourth sensor wasn't efficient!
                                //but only if it could have been
                                if (!s4->isDeadArea(t->findXYPrediction(s4->zPosition()))) {
                                        s4->registerTrackMiss(t->fourthSensorThresh(),
                                                        t->findXYPrediction(s4->zPosition()),
                                                        t->timeStamp());
                                        ineff++;
                                        //Where would the hit have been i.t.o. a pixel coordinate?
                                        MapsSensor::coord predHit(0, 0);
                                        s4->convertPhysicalToHit(
                                                        t->findXYPrediction(s4->zPosition()), predHit);
                                        if (predHit.first < 84)
                                                shapIneff++;
                                        else
                                                sampIneff++;
                                } else {
                                        //hit would have intersected with a dead area
                                        ++deadArs;
                                }
                        }

                }
        }

        //All that follows is housekeeping and plot extraction...

        char output[100];
        strcpy(output, argv[2]);
        TFile f(output, "recreate");

        //Find the projected dead area defined by the 4 sensors
        TH2F p;
        MapsSensor::mutualSensorMask(sensors, p);
        
        f.mkdir("effCurves");
        
        //iterate over sensors to extract efficiency plots
        for (std::vector<MapsSensor*>::iterator it = sensors.begin(); it
                        != sensors.end(); it++) {
                MapsSensor* s = *it;
                TH1F h, r, q;
                TH1F g, k;
                TH2F j, l, m, n;
                f.cd("effCurves");
                //Get global sensor efficiency
                s->getEfficiencyCurve(h, 13, 80, 210, true, true);
                h.Write();
                //Just shapers please
                s->getEfficiencyCurve(r, 13, 80, 210, true, false);
                r.Write();
                //Now just samplers
                s->getEfficiencyCurve(q, 13, 80, 210, false, true);
                q.Write();
                f.cd("/");

                s->getEfficiencyTimestamps(g);
                s->getEfficiencyXYPlot(j);
                s->getInefficiencyXYPlot(m, true);
                s->getInefficiencyTimestamps(k);
                s->getFourthHitResidualPlot(l);
                s->getEfficiencyByGroup(n);

                g.Write();
                j.Write();
                k.Write();
                l.Write();
                m.Write();
                n.Write();
                p.Write();
                std::cout.width(4);
                std::cout << "\n" << *s << ": Efficiency: \n";
                //              for (unsigned tMult(7); tMult < 20; tMult++) {
                //                      double efficiency = s->getEfficiency(tMult*10);
                //                      double shapEff = s->getEfficiency(tMult*10, true, false);
                //                      double sampEff = s->getEfficiency(tMult*10, false, true);
                //                      if (efficiency > -0.9) {
                //                              std::cout.precision(3);
                //
                //                              std::cout << "\t Thresh: " << tMult*10 << "\t: " << efficiency
                //                                              * 100 << ",\t shap:\t" << shapEff * 100
                //                                              << ",\t samp:\t" << sampEff * 100 << "\n";
                //                      }
                //              }
                
                //Commented code above now obseleted by...
                s->printEfficiencies(std::cout);
        }

        //save plots to disk
        std::cout << "\tWriting root file: " << output << "\n";
        f.Write();
        f.Close();

        //close up shop
        std::cout << "ExtractEfficiencies: summary:\n";
        std::cout << "\tTotal candidate tracks:\t" << candidateTracks << "\n";
        std::cout << "\tEfficient hits:\t" << eff << "\n";
        std::cout << "\tInefficient hits:\t" << ineff << "\n";
        std::cout << "\t\tShaper efficient hits:\t" << shapEff << "\n";
        std::cout << "\t\tShaper inefficient hits:\t" << shapIneff << "\n";
        std::cout << "\t\tSampler efficient hits:\t" << sampEff << "\n";
        std::cout << "\t\tSampler inefficient hits:\t" << sampIneff << "\n";
        std::cout << "\tDead area intersections:\t" << deadArs << "\n";
        std::cout << "Global efficiency:\t" << 100 * static_cast<double>(eff)
                        / (eff+ineff) << "\n";
        std::cout << "\tShapers:\t" << 100 * static_cast<double>(shapEff)
                        / (shapEff+shapIneff) << "\n";
        std::cout << "\tSamplers:\t" << 100 * static_cast<double>(sampEff)
                        / (sampEff+sampIneff) << "\n";

        std::cout << "Have a nice day.\n";

}

---