FigureOfMerit.cxx Source File

00001 //  $Header: /nfs/slac/g/glast/ground/cvs/merit/src/FigureOfMerit.cxx,v 1.7 2001/11/01 17:28:00 burnett Exp $
00002 
00003 
00004 #include "FigureOfMerit.h"
00005 #include "Cut.h"
00006 #include "MultiPSF.h"
00007 #include "Level1.h"
00008 
00009 #include <cmath>
00010 #include <iomanip>
00011 #include <algorithm>
00012 
00013 // put these here because gcc could not find the inline versions???
00014 unsigned        FigureOfMerit::accepted() const { return m_accepted; }
00015 float    FigureOfMerit::area(){return s_area;}
00016 unsigned FigureOfMerit::generated()  { return s_generated; }
00017 
00018 using namespace std;
00019 
00020 static inline double sqr(double x){return x*x;}
00021 //=============================================================================
00022 // Figure of Merit analysis cut functions: allow dynamic allocation/ordering of a list
00023 //      of analysis cuts which directly interface with Figure of Merit.
00024 
00025 //=============================================================================
00026 // utility that sums the variable
00027 class Summation : public Analyze {
00028 public:
00029     Summation(const Tuple&t , const char* item, const char* label )
00030         : Analyze(t, item, label), m_total(0) {};
00031     float average()const{return count() ? m_total/count():0;}
00032     float total()const{return m_total;}
00033 private:
00034     float  m_total; 
00035     virtual bool apply () { m_total += item();  return true;   }
00036 };
00037 //=============================================================================
00038 class FOMelapsed : public Analyze {
00039 public:
00040     FOMelapsed(const Tuple&t ): Analyze(t, "Triage_Time", "Elapsed time (sec):"), m_total(0),m_last(0) {};
00041     void FOMelapsed::report(ostream& out)
00042     {
00043         out << endl << make_label(name());
00044         out << setw(6) << m_total;
00045     }
00046     double time()const{return m_total;}
00047 
00048 private:
00049     float  m_total, m_last;
00050 
00051     virtual bool apply (){ 
00052         float val = item();
00053         if( val < m_last ) m_last=0;  // in case concatenated input files
00054         m_total += val - m_last;
00055         m_last = val;
00056         return    true;
00057     };
00058 };
00059 
00060 //=============================================================================
00061 class FOMtrigrate : public FOMelapsed {
00062 public:
00063     FOMtrigrate(const Tuple&t ): FOMelapsed(t) {};
00064     void FOMtrigrate::report(ostream& out)
00065     {
00066         out << endl << make_label("Rate:");
00067         out << setw(6) << setprecision(3);
00068         float rate = (time() ? count()/time() : 0);
00069         if( rate<1000. )  out  << rate << " Hz";
00070         else out << rate/1000. << " kHz";
00071     }
00072 
00073 private:
00074 };
00075 
00076 //=============================================================================
00077 class FOMdeadtime : public Summation {
00078 public:
00079     FOMdeadtime(const Tuple&t ): Summation(t, "Dead_Time", "Dead Time:"){};
00080     void FOMdeadtime::report(ostream& out)
00081     {
00082         out << endl << make_label(name());
00083         out << setw(6) << total();
00084     }
00085 };
00086 //=============================================================================
00087 class FOMROtime : public Summation {
00088 public:
00089     FOMROtime(const Tuple&t ): Summation(t, "Max_RO_Time", "Readout Time (max):"){};
00090     void FOMROtime::report(ostream& out)
00091     {
00092         out << endl << make_label(name());
00093         out << setw(6) << average();
00094     }
00095 };
00096 
00097 //=============================================================================
00098 class FOMevtsize : public Summation {
00099 public:
00100     FOMevtsize(const Tuple&t ): Summation(t, "Total_Evt_Size", "Event Size (total):"){};
00101     void FOMevtsize::report(ostream& out)
00102     {
00103         out << endl << make_label(name());
00104         out << setw(6) << average() << " bits";
00105     }
00106 };
00107   //=============================================================================
00108 class EventSize : public Analyze{
00109 public:
00110     EventSize(const Tuple&t ) : Analyze("Average Event size"),
00111         m_acd(t, "ACD_TileCount", "ACD Hits (19 b)"),
00112         m_ssd(t, "TKR_Cnv_Lyr_Hits", "SSD Hits (20 b)"),
00113         m_cal(t, "Cal_No_Xtals", "CAL Hits (40 b)")
00114     {}
00115     void EventSize::report(ostream& out)
00116     {   
00117         float size 
00118             = 19* m_acd.stat().mean()
00119             + 20* m_ssd.stat().mean()
00120             + 40* m_cal.stat().mean();
00121         separator(out);
00122         out << endl << make_label(name()) 
00123             << setw(6)<< setprecision(3) << size/1e3 <<" kbits";
00124 
00125         m_acd.report(out);
00126         m_ssd.report(out);
00127         m_cal.report(out);
00128     }
00129     bool apply(){
00130         m_acd(); m_ssd(); m_cal();
00131     return true;}
00132 
00133 private:
00134     Statistic m_acd;
00135     Statistic m_ssd;
00136     Statistic m_cal; 
00137 };
00138 //=============================================================================
00139 class FOMnsistrips : public Summation {
00140 public:
00141     FOMnsistrips(const Tuple&t ): Summation(t, "TKR_Cnv_Lyr_Hits", "Tracker hits (avg):"){};
00142     void FOMnsistrips::report(ostream& out)
00143     {
00144         out << endl << make_label(name());
00145         out << setw(6) << average() << " hits";
00146     }
00147 
00148 };
00149 
00150 //=============================================================================
00151 class FOMncsixtals : public Summation {
00152 public:
00153     FOMncsixtals(const Tuple&t ): Summation(t, "Cal_No_Xtals", "CsI logs hit (avg):"){};
00154     void FOMncsixtals::report(ostream& out)
00155     {
00156         out << endl << make_label(name());
00157         out << setw(6) << average() << " logs";
00158     }
00159 };
00160 
00161 //=============================================================================
00162 class FOMnacdtiles : public Summation {
00163 public:
00164     FOMnacdtiles(const Tuple&t ): Summation(t, "ACD_TileCount", "ACD tiles hit (avg):"){};
00165     void FOMnacdtiles::report(ostream& out)
00166     {
00167         out << endl << make_label(name());
00168         out << setw(6) << average() << " tiles";
00169     }
00170 };
00171 
00172 //=============================================================================
00173 class FOMtrig1Cal : public Analyze {
00174 public:
00175     FOMtrig1Cal(const Tuple&t, int mask=64 )
00176         : Analyze(t, "Trig_Bits", "Level1")
00177         ,m_mask(mask){};
00178 private:
00179     virtual bool apply () {
00180         return ( (static_cast<int>(item()) & m_mask )==m_mask); } //  cal
00181     int m_mask;
00182 };
00183 
00184 //=============================================================================
00185 // Level  2 analysis
00186 class Level2 : public Analyze { 
00187 public:
00188     Level2(const Tuple& t): Analyze(t, "Trig_Level2_Bits", " Level2")
00189         , m_trig_cal(t,32)
00190     { clear(); };
00191     void report(ostream& out)
00192     {
00193         out << endl << name() ;
00194         out << endl << "          track:" << setw(8) << m_track;
00195         out << endl << "           veto:" << setw(8) << m_veto;
00196         out << endl << "  veto & !hical:" << setw(8) << m_cal;
00197         Analyze::report(out);
00198         separator(out);
00199     }
00200     void clear(){m_track=m_veto=m_cal=0; Analyze::clear();}
00201 private:
00202     virtual bool  apply ()
00203     { 
00204         int trig2 = static_cast<unsigned>(item()) & 15;
00205         bool cal = m_trig_cal() ;
00206         if( trig2 & 1){
00207             m_track++; 
00208             if( trig2 >1 ){
00209                 m_veto++;
00210                 if( !cal) m_cal++ ;
00211             }
00212         }
00213         return trig2 == 1 ||  cal;
00214     }
00215     FOMtrig1Cal m_trig_cal;
00216     int m_track; // number with track bit set
00217     int m_veto;  // number with also veto
00218     int m_cal;   // calorimeter trigger overrides veto
00219 
00220 };
00221 //=============================================================================
00222 class FOMtrigII : public Analyze {      
00223 public:
00224     FOMtrigII(const Tuple& t): Analyze(t, "Trig_Level2_Bits", "Level2 3row only"){};
00225 private:
00226     virtual bool  apply (){ return ((static_cast<unsigned>(item()) & 1) == 1); }
00227 };
00228 //=============================================================================
00229 class FOML2T : public Analyze {
00230 public:
00231     FOML2T (const Tuple& t) 
00232         : Analyze (t, "ACD_DOCA", "L2 Selections")
00233         , m_sourceid(t.tupleItem ("MC_src_Id"))
00234           , m_numtracks(t.tupleItem ("TKR_No_Tracks"))
00235           , m_edeposit(t.tupleItem ("Cal_Energy_Deposit"))
00236 
00237     {    }
00238 private:
00239     virtual bool    apply () {
00240         float doca = item();
00241             float sourceid = *m_sourceid;
00242         float numtracks = *m_numtracks;
00243         float edeposit = *m_edeposit;
00244 
00245         return (sourceid!=2&&((numtracks>0.&&doca>25.)||(edeposit>10000.)));
00246     }
00247     
00248     const TupleItem*    m_sourceid;
00249     const TupleItem*    m_numtracks;
00250     const TupleItem*    m_edeposit;
00251 };
00252 //=============================================================================
00253 class FOML3T : public Analyze {
00254 public:
00255     FOML3T (const Tuple& t) 
00256         : Analyze (t, "Cal_Energy_Deposit", "L3 Selections")
00257         , m_surplushit(t.tupleItem ("REC_Surplus_Hit_Ratio"))
00258           , m_xtalrat(t.tupleItem ("Cal_Xtal_Ratio"))
00259           , m_fiterrnrm(t.tupleItem ("Cal_Fit_errNrm"))
00260 
00261     {    }
00262 private:
00263     virtual bool    apply () {
00264         float edeposit = item();
00265         float surplushit = *m_surplushit;
00266         float xtalrat = *m_xtalrat;
00267         float fiterrnrm = *m_fiterrnrm;
00268 
00269         return (((edeposit<1.&&surplushit>2.)||(xtalrat>.2))&&fiterrnrm<15.); 
00270 //        return (xtalrat>.2&&fiterrnrm<15.); 
00271     }
00272     const TupleItem*    m_surplushit;
00273     const TupleItem*    m_xtalrat;
00274     const TupleItem*    m_fiterrnrm;
00275 };
00276 //=============================================================================
00277 class FOML1V : public Analyze {
00278 public:
00279     FOML1V (const Tuple& t) 
00280         : Analyze (t, "ACD_Throttle_Bits", "L1 Throttle")
00281         , m_tilecount(t.tupleItem ("ACD_TileCount"))
00282           , m_siderow3(t.tupleItem ("ACD_No_SideRow3"))
00283           , m_siderow2(t.tupleItem ("ACD_No_SideRow2"))
00284           , m_trigbits(t.tupleItem ("Trig_Bits"))
00285 
00286     {    }
00287 private:
00288     virtual bool    apply () {
00289         float throttlebits = item();
00290           float tilecount = *m_tilecount;
00291         float siderow3 = *m_siderow3;
00292         float siderow2 = *m_siderow2;
00293         int trigbits = *m_trigbits;
00294 
00295         return (((throttlebits==0||throttlebits>127.)&&(tilecount-siderow2-siderow3)<3.)||(trigbits&16));
00296     }
00297     
00298     const TupleItem*    m_throttlebits;
00299     const TupleItem*    m_tilecount;
00300     const TupleItem*    m_siderow3;
00301     const TupleItem*    m_siderow2;
00302     const TupleItem*    m_trigbits;
00303 };
00304 //=============================================================================
00305 class WriteTuple : public Analyze {
00306 public:
00307     WriteTuple(const Tuple& t) : Analyze( "Written to output tuple"),m_tuple(t)
00308     {
00309         m_tuple.writeHeader(std::cout);
00310     }
00311 private:
00312     virtual bool    apply () {std::cout << m_tuple; ; return true; }
00313     const Tuple& m_tuple;
00314 };
00315 //=============================================================================
00316 // Only the Front
00317 class FOMFront : public AnalysisList {  
00318 public: 
00319 
00320     FOMFront(const Tuple& t): AnalysisList(" TKR Front only")
00321     {
00322         push_back( new Cut(t, "TKR_First_XHit<12.") );
00323     };
00324     void report(ostream& out)
00325     {
00326         separator(out);
00327         AnalysisList::report(out);
00328     }
00329 };
00330 //=============================================================================
00331 // Only the Back
00332 class FOMBack : public AnalysisList {   
00333 public: 
00334 
00335     FOMBack(const Tuple& t): AnalysisList(" TKR Back only")
00336     {
00337         push_back( new Cut(t, "TKR_First_XHit>11.") );
00338     };
00339     void report(ostream& out)
00340     {
00341         separator(out);
00342         AnalysisList::report(out);
00343     }
00344 };
00345 //=============================================================================
00346 // Analyze the level 3 cuts
00347 class Level3 : public AnalysisList {    
00348 public: 
00349 
00350     Level3(const Tuple& t): AnalysisList(" Level 3")
00351     {
00352         push_back( new Cut(t, "Cal_No_Xtals>0") );
00353         push_back( new Cut(t, "TKR_No_Tracks>0") );
00354         push_back( new Cut(t, "ACD_DOCA>25") );
00355         //push_back( new Cut(t, "CsI_Fit_errNrm>10"));
00356     };
00357     void report(ostream& out)
00358     {
00359         separator(out);
00360         AnalysisList::report(out);
00361     }
00362 
00363 };
00364 //=============================================================================
00365 // Analyze the Atwood cuts
00366 class CosmicCuts : public AnalysisList {
00367 public:
00368     CosmicCuts(const Tuple&t, bool noline=false)
00369         : AnalysisList("  ---Cosmic cuts---", noline)
00370     {
00371         push_back( new Cut(t, "REC_Surplus_Hit_Ratio", Cut::GT, 2.05) );
00372         push_back( new Cut(t, "ACD_DOCA",      Cut::GT, 30) );
00373         push_back( new Cut(t, "Cal_Fit_errNrm", Cut::LT, 5) );
00374         push_back( new Cut(t, "Cal_Xtal_Ratio",  Cut::GT, 0.25) );
00375     }
00376 };
00377 
00378 
00379 //=============================================================================
00380 class FOMaccepted : public Analyze {
00381 public:
00382     FOMaccepted(FigureOfMerit* f) : Analyze("Accepted for analysis"),m_f(f) {}
00383 private:
00384     virtual bool    apply () { m_f->accept(); return true; }
00385     FigureOfMerit* m_f;
00386 };
00387 
00388 //=============================================================================
00389 
00390 const Tuple*     FigureOfMerit::s_tuple=0;
00391 
00392 
00393 unsigned FigureOfMerit::s_generated = 100000;
00394 double FigureOfMerit::s_area = 60000.;
00395 
00396 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00397 void    FigureOfMerit::setCuts ( string istr )
00398 {
00399     delete m_cuts;
00400     m_cuts = new AnalysisList( string("Analysis cuts: ")+istr );
00401     string::const_iterator      it = istr.begin();
00402     
00403     m_cuts->push_back( new Analyze("Found in tuple") );  // an event count for starting (# in tuple)
00404     m_cuts->push_back( new Statistic(*s_tuple, "MC_Energy", "Generated energy"));
00405     
00406     while (it != istr.end()) {
00407         switch (*it++) {
00408             
00409         case 'G': // Gnn, -- override number generated
00410             s_generated= ::atoi(it);
00411             while( it !=istr.end() && (*it++)!=',');
00412             cerr << "Overriding generated number to " << s_generated << endl;
00413 
00414             break;
00415         case '1':           // 1(one) = level I trigger
00416             m_cuts->push_back( new Level1(*s_tuple) );
00417             break;
00418             // added 18 oct 2001 by S.Ritz
00419         case 'V':        // V = level 1 VETO using ACD
00420             m_cuts->push_back( new FOML1V(*s_tuple) );
00421             break;
00422             
00423             //change 18 oct 2001 by S.Ritz to new L2T and L3T selections
00424         case '2':           // 2 = level II trigger - no veto
00425             m_cuts->push_back( new FOML2T(*s_tuple) );
00426             break;
00427         case '3':           // 3 = level 3 trigger 
00428             m_cuts->push_back( new FOML3T(*s_tuple) );
00429             break;
00430             // 18 october 2001 by S.Ritz -- add selections for Front and Back
00431         case 'F':               // F = Front TKR section only
00432             m_cuts->push_back( new FOMFront(*s_tuple) );
00433             break;
00434         case 'B':               // B = Back TKR section only
00435             m_cuts->push_back( new FOMBack(*s_tuple) );
00436             break;
00437         case 'd':               // d = level II trigger - with veto
00438             m_cuts->push_back( new FOMtrigII(*s_tuple) );
00439             break;
00440         case 'n':           // n = ntracks
00441             m_cuts->push_back( new Cut(*s_tuple, "TKR_No_Tracks>0" ) );
00442             break;
00443         case 'c':           // c = Cosmic cuts
00444             m_cuts->push_back( new CosmicCuts(*s_tuple) );
00445             break;
00446         case 's':
00447             m_cuts->push_back( new EventSize(*s_tuple) );
00448             break;
00449         case 't':
00450             m_cuts->push_back( new FOMROtime(*s_tuple) );
00451             break;
00452         case 'M': //M0, M1, ...
00453             m_cuts->push_back( new MultiPSF(*s_tuple, *(it++)) );
00454             break;
00455             
00456         case 'P': /* P = PSF analysis   */  m_cuts->push_back( new PSFanalysis(*s_tuple) );   break;
00457         case 'E':           break;
00458             
00459         case 'A': /* A = accepted */
00460             m_layers.push_back(LayerGroup(*s_tuple,0,11));
00461             m_layers.push_back(LayerGroup(*s_tuple,12,15));
00462             
00463             m_cuts->push_back( new FOMaccepted(this) );     break;
00464             
00465         case 'W': /* W = Write */        m_cuts->push_back( new WriteTuple(*s_tuple)); break;
00466         case 'D': /* D = dead time */    m_cuts->push_back( new FOMdeadtime(*s_tuple) ); break;
00467         case 'L': /* L = elapsed time */ m_cuts->push_back( new FOMelapsed(*s_tuple) ); break;
00468         case 'R': /* R = trigger rate */ m_cuts->push_back( new FOMtrigrate(*s_tuple) ); break;
00469         case '(': // (cut) -- pass in the iterator
00470             m_cuts->push_back(new Cut(*s_tuple, it, istr.end() )); break; 
00471         case 'X': // Xname, statistic on name
00472             m_cuts->push_back(new Statistic(*s_tuple, it, istr.end() )); break;
00473         default: 
00474             cerr << "Key '" << *(it-1) << "' ignored" << endl;
00475               break;
00476         }   // switch
00477     }   // while
00478     
00479 }
00480 //=============================================================================
00481 FigureOfMerit::FigureOfMerit(const Tuple& t, std::string cut_string)
00482 : m_cuts(0)
00483 , m_accepted(0)
00484 {
00485     if (!s_tuple)       {   // initialize static member functions
00486         s_tuple = &t;
00487 
00488         // analyze title for parameters
00489         string title(t.title());
00490         string::size_type pos = title.find("area =");
00491         if( pos != string::npos ) {
00492             string a(title, pos+6, title.length() );
00493             s_area = atof(a.data());
00494         }
00495         else {
00496             cerr << "Area not found in title: assuming "
00497                 << s_area <<" cm^2" << endl;
00498         }
00499 
00500         // look for number generated in title: either generated:1000 or gen(1000)
00501         pos = title.find("generated");
00502         if( pos != string::npos ) {
00503             string a(title, pos+10, 10);
00504             s_generated = atoi(a.data());
00505         }
00506         else if ( (pos = title.find("gen(")) != string::npos) {
00507             string a(title, pos+4, 10);
00508             s_generated = atoi(a.data());
00509         }
00510         else    cerr << "generated events not found in title" << endl;
00511     }   // if (!s_instance)
00512 
00513     if( !cut_string.empty()) setCuts(cut_string);
00514 }
00515 
00516 
00517 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00518 void FigureOfMerit::accept()
00519 {
00520     m_accepted++;
00521     for( std::vector<LayerGroup>::iterator layer = m_layers.begin(); layer!=m_layers.end(); ++layer) (*layer)();
00522 
00523 }
00524 
00525 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00526 void FigureOfMerit::execute()
00527 {
00528     (*m_cuts)();
00529 }
00530 
00531 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00532 void FigureOfMerit::report(ostream & out)
00533 {
00534     out << "\n"<< Analyze::make_label("Generated events") << setw(6) << generated() ;
00535 
00536     Analyze::separator(out,'=');
00537     m_cuts->report(out);        // print analysis cuts
00538 
00539     if( ! accepted() ) return;
00540 
00541     float  area_per_event = s_area/s_generated,  fom_total=0;
00542 
00543     // loop over the ranges of layers
00544     for( std::vector<LayerGroup>::iterator layers = m_layers.begin(); layers!=m_layers.end(); ++layers){
00545         layers->report(out);
00546         float eff_area = area_per_event*layers->count();
00547         
00548         out     << "\n" << Analyze::make_label(" effective area")
00549             << setw(6)
00550             << static_cast<int>(eff_area+0.5) << " cm^2";
00551         float fom = sqrt(eff_area)/layers->sigma();
00552         out << "\n" << Analyze::make_label("Figure of merit")
00553             << setw(6)
00554             << static_cast<int>(fom+0.5) << " cm";
00555         
00556         Analyze::separator(out);
00557         fom_total = sqrt(sqr(fom_total)+sqr(fom));
00558     }
00559     out << "\n" << Analyze::make_label("total effective area")
00560         << setw(6)
00561         << static_cast<int>(area_per_event*accepted()+0.5) << " cm^2";
00562 
00563     out << "\n" << Analyze::make_label("Combined FOM")
00564         << setw(6) << static_cast<int>(fom_total+0.5) << " cm";
00565 
00566     out << endl;
00567 }
00568 
00569 
00570