tests/src/ExtractEfficienciesCosmics.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 <iostream>
#include "TRandom2.h"
#include "TFile.h"
#include "TH1F.h"
#include "Operators.h"

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.1, 0.1);
        //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();

        //count efficient and inefficient frequency
        unsigned ineff(0), eff(0), deadArs(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;
                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++;
                        } else {
                                //fourth sensor wasn't efficien!
                                //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++;
                                } else {
                                        ++deadArs;
                                }
                        }

                }
        }
        
        //All that follows is housekeeping and plot extraction...
        
        char output[100];
        strcpy(output, argv[2]);
        TFile f(output, "recreate");

        //iterate over sensors to extract efficiency plots      
        std::vector<MapsSensor*>& sensors = mtm2.getSensors();
        for (std::vector<MapsSensor*>::iterator it = sensors.begin(); it
                        != sensors.end(); it++) {
                MapsSensor* s = *it;
                TH1F h;
                TH1F g, k;
                TH2F j, l, m, n;
                s->getEfficiencyCurve(h, 20, 0, 200);
                s->getEfficiencyTimestamps(g);
                s->getEfficiencyXYPlot(j);
                s->getInefficiencyXYPlot(m, true);
                s->getInefficiencyTimestamps(k);
                s->getFourthHitResidualPlot(l);
                s->getEfficiencyByGroup(n);
                h.Write();
                g.Write();
                j.Write();
                k.Write();
                l.Write();
                m.Write();
                n.Write();
                std::cout.precision(4);
                std::cout << "\n" << *s << ": Efficiency: \n";
                for (unsigned tMult(7); tMult < 20; tMult++) {
                        double efficiency = s->getEfficiency(tMult*10);
                        if (efficiency > -0.9)
                                std::cout << "\t Thresh: " << tMult*10 << "\t: " << efficiency
                                                * 100 << "\n";
                }
        }

        //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 << "\tDead area intersections:\t" << deadArs << "\n";

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

}

---