tests/src/ExtractEfficiencies.cpp

Go to the documentation of this file.
00001 #include "MapsTrackManager.hh"
00002 #include "MapsTrack.hh"
00003 #include "MapsSensor.hh"
00004 #include "MapsException.hh"
00005 #include <cmath>
00006 #include <algorithm>
00007 #include <functional>
00008 #include <iostream>
00009 #include "TRandom2.h"
00010 #include "TFile.h"
00011 #include "TH1F.h"
00012 #include "Operators.h"
00013 
00014 /** ExtractEfficiencies.cpp
00015  * 
00016  * Workhorse for determining sensor efficiency once tracks have been processed
00017  * and you're happy with the alignment.
00018  * 
00019  * Exports plots to supplied outpufile name.
00020  * 
00021  * This is a good place to start to understand the general code structure, and
00022  * as a general example for how you could write your own tool.
00023  * 
00024  * Usage: ExtractEfficiencies input file outputfile
00025  * 
00026  * Jamie Ballin, Imperial College London
00027  *              February 2008
00028  * 
00029  */
00030 
00031 int main(int argc, const char **argv) {
00032 
00033         if (argc != 3) {
00034                 std::cout << "Usage: ExtractEfficiencies <input file> <outputfile>\n";
00035                 return 0;
00036         }
00037         std::cout << "Welcome to ExtractEfficiencies...\n";
00038 
00039         //Define chi squared probability cuts
00040         double cut(0.05);
00041         //Define errors for chi squared calculation
00042         std::pair<double, double> error(0.018, 0.018);
00043         //std::pair<double, double> error(0.026, 0.02);
00044         std::cout << "\tUsing cuts for probability of: \t" << cut << "\n";
00045         std::cout << "\tError parameters: \t" << error << "\n";
00046 
00047         std::cout << "\tRecreating from root file...\n";
00048         char input[100];
00049         strcpy(input, argv[1]);
00050 
00051         //Create a maps track manager to recreate the tracks
00052         MapsTrackManager mtm2;
00053         //Reinitialise all the tracks and sensors in mtm2
00054         mtm2.recreateFromRootFile(input);
00055 
00056         std::vector<MapsTrack*>& tracks = mtm2.getTracks();
00057         std::vector<MapsSensor*>& sensors = mtm2.getSensors();
00058         MapsSensor* s1;
00059         MapsSensor* s2;
00060         for (std::vector<MapsSensor*>::const_iterator it = sensors.begin(); it
00061                         != sensors.end(); ++it) {
00062                 if ((*it)->id() == 8) {
00063                         s1 = *it;
00064                 } else if ((*it)->id() == 6) {
00065                         s2 = *it;
00066                 }
00067         }
00068 
00069         //count efficient and inefficient frequency
00070         unsigned ineff(0), eff(0), deadArs(0);
00071         //counts specific to shapers and samplers
00072         unsigned shapEff(0), sampEff(0);
00073         unsigned shapIneff(0), sampIneff(0);
00074         unsigned candidateTracks(0);
00075 
00076         for (std::vector<MapsTrack*>::iterator it = tracks.begin(); it
00077                         != tracks.end(); it++) {
00078                 MapsTrack* t = *it;
00079 
00080                 //this is not obligatory, but highly advised!
00081                 t->setTrackError(error);
00082 
00083                 //see if we have a good three hit track first: for tracks which are just three hits anyway
00084                 //t3 will be a copy of t
00085                 MapsTrack t3;
00086                 //construct the three hit track
00087                 t->make3HitTrack(t3);
00088                 t3.setTrackError(error);
00089 
00090                 //evaluate quality of the 3 hit track
00091                 if (t3.chiSqProb(0) > cut && t3.chiSqProb(1) > cut) {
00092                         //it's a real track then
00093                         candidateTracks++;
00094 
00095                         //get the original track's fourth sensor
00096                         MapsSensor* s4 = t->fourthSensor();
00097 
00098                         //pose the ultimate question...
00099                         std::pair<double, double> resid;
00100                         unsigned fourthConfirm = t->getFourthHitResidual(t3, resid);
00101 
00102                         if (fourthConfirm == 0 && t->chiSqProb(0) > cut && t->chiSqProb(1)
00103                                         > cut) {
00104                                 //fourth sensor was efficient
00105                                 //one could also cut on the residual value, if you wanted
00106                                 s4->registerTrackConfirm(t->fourthSensorThresh(), t->getHits().at(s4), t->timeStamp(), resid);
00107                                 eff++;
00108                                 if (t->getHits().at(s4).first < 84)
00109                                         shapEff++;
00110                                 else
00111                                         sampEff++;
00112                         } else {
00113                                 //fourth sensor wasn't efficient!
00114                                 //but only if it could have been
00115                                 if (!s4->isDeadArea(t->findXYPrediction(s4->zPosition()))) {
00116                                         s4->registerTrackMiss(t->fourthSensorThresh(),
00117                                                         t->findXYPrediction(s4->zPosition()),
00118                                                         t->timeStamp());
00119                                         ineff++;
00120                                         //Where would the hit have been i.t.o. a pixel coordinate?
00121                                         MapsSensor::coord predHit(0, 0);
00122                                         s4->convertPhysicalToHit(
00123                                                         t->findXYPrediction(s4->zPosition()), predHit);
00124                                         if (predHit.first < 84)
00125                                                 shapIneff++;
00126                                         else
00127                                                 sampIneff++;
00128                                 } else {
00129                                         //hit would have intersected with a dead area
00130                                         ++deadArs;
00131                                 }
00132                         }
00133 
00134                 }
00135         }
00136 
00137         //All that follows is housekeeping and plot extraction...
00138 
00139         char output[100];
00140         strcpy(output, argv[2]);
00141         TFile f(output, "recreate");
00142 
00143         //Find the projected dead area defined by the 4 sensors
00144         TH2F p;
00145         MapsSensor::mutualSensorMask(sensors, p);
00146         
00147         f.mkdir("effCurves");
00148         
00149         //iterate over sensors to extract efficiency plots
00150         for (std::vector<MapsSensor*>::iterator it = sensors.begin(); it
00151                         != sensors.end(); it++) {
00152                 MapsSensor* s = *it;
00153                 TH1F h, r, q;
00154                 TH1F g, k;
00155                 TH2F j, l, m, n;
00156                 f.cd("effCurves");
00157                 //Get global sensor efficiency
00158                 s->getEfficiencyCurve(h, 13, 80, 210, true, true);
00159                 h.Write();
00160                 //Just shapers please
00161                 s->getEfficiencyCurve(r, 13, 80, 210, true, false);
00162                 r.Write();
00163                 //Now just samplers
00164                 s->getEfficiencyCurve(q, 13, 80, 210, false, true);
00165                 q.Write();
00166                 f.cd("/");
00167 
00168                 s->getEfficiencyTimestamps(g);
00169                 s->getEfficiencyXYPlot(j);
00170                 s->getInefficiencyXYPlot(m, true);
00171                 s->getInefficiencyTimestamps(k);
00172                 s->getFourthHitResidualPlot(l);
00173                 s->getEfficiencyByGroup(n);
00174 
00175                 g.Write();
00176                 j.Write();
00177                 k.Write();
00178                 l.Write();
00179                 m.Write();
00180                 n.Write();
00181                 p.Write();
00182                 std::cout.width(4);
00183                 std::cout << "\n" << *s << ": Efficiency: \n";
00184                 //              for (unsigned tMult(7); tMult < 20; tMult++) {
00185                 //                      double efficiency = s->getEfficiency(tMult*10);
00186                 //                      double shapEff = s->getEfficiency(tMult*10, true, false);
00187                 //                      double sampEff = s->getEfficiency(tMult*10, false, true);
00188                 //                      if (efficiency > -0.9) {
00189                 //                              std::cout.precision(3);
00190                 //
00191                 //                              std::cout << "\t Thresh: " << tMult*10 << "\t: " << efficiency
00192                 //                                              * 100 << ",\t shap:\t" << shapEff * 100
00193                 //                                              << ",\t samp:\t" << sampEff * 100 << "\n";
00194                 //                      }
00195                 //              }
00196                 
00197                 //Commented code above now obseleted by...
00198                 s->printEfficiencies(std::cout);
00199         }
00200 
00201         //save plots to disk
00202         std::cout << "\tWriting root file: " << output << "\n";
00203         f.Write();
00204         f.Close();
00205 
00206         //close up shop
00207         std::cout << "ExtractEfficiencies: summary:\n";
00208         std::cout << "\tTotal candidate tracks:\t" << candidateTracks << "\n";
00209         std::cout << "\tEfficient hits:\t" << eff << "\n";
00210         std::cout << "\tInefficient hits:\t" << ineff << "\n";
00211         std::cout << "\t\tShaper efficient hits:\t" << shapEff << "\n";
00212         std::cout << "\t\tShaper inefficient hits:\t" << shapIneff << "\n";
00213         std::cout << "\t\tSampler efficient hits:\t" << sampEff << "\n";
00214         std::cout << "\t\tSampler inefficient hits:\t" << sampIneff << "\n";
00215         std::cout << "\tDead area intersections:\t" << deadArs << "\n";
00216         std::cout << "Global efficiency:\t" << 100 * static_cast<double>(eff)
00217                         / (eff+ineff) << "\n";
00218         std::cout << "\tShapers:\t" << 100 * static_cast<double>(shapEff)
00219                         / (shapEff+shapIneff) << "\n";
00220         std::cout << "\tSamplers:\t" << 100 * static_cast<double>(sampEff)
00221                         / (sampEff+sampIneff) << "\n";
00222 
00223         std::cout << "Have a nice day.\n";
00224 
00225 }
00226 

Generated on Wed Mar 19 17:47:58 2008 for MapsTracks by  doxygen 1.5.2