GlastRecon.cxx Source File

00001 //   $Header: /nfs/slac/g/glast/ground/cvs/reconstruction/src/GlastRecon.cxx,v 1.11 2000/12/19 03:30:14 igable Exp $  
00002 
00003 #ifdef __GNUG__                                    
00004 #pragma implementation
00005 #endif
00006 
00007 #include "reconstruction/GlastRecon.h"
00008 
00009 #include "reconstruction/MCRecon.h"
00010 #include "reconstruction/CalRecon.h"
00011 #include "reconstruction/TrackerRecon.h"
00012 #include "reconstruction/TriggerRecon.h"
00013 #include "reconstruction/VetoRecon.h"
00014 #include "reconstruction/CalProfile.h"
00015 #include "reconstruction/ReconVisitor.h"
00016 
00017 #include "DetectorData.h"
00018 
00019 //THB #include "graphics/PrintControl.h"
00020 
00021 #include "data/GlastData.h"
00022 
00023 
00024 GlastRecon::GlastRecon ()
00025 :Recon()
00026 {
00027     setTitle("GlastRecon: ver 2.0");
00028 //TODO: Remove this hack
00029 
00030    
00031     // Note: the order of reconstruction depends on the following ordering
00032     addRecon(mc   = new MCRecon());
00033     addRecon(trig = new TriggerRecon());
00034     addRecon(cal  = new CalRecon);
00035     addRecon(trkr = new TrackerRecon(this));
00036     addRecon(veto = new VetoRecon( cal, trkr ));
00037     addRecon(prof = new CalProfile( cal, trkr ));
00038     
00039     i_gl_xfit_err   = data.addElement("Fit_Xdir_Err");
00040     i_gl_yfit_err   = data.addElement("Fit_Ydir_Err");
00041     i_gl_zfit_err   = data.addElement("Fit_Zdir_Err");
00042     i_gl_xGam_err   = data.addElement("Gamma_Xdir_Err");
00043     i_gl_yGam_err   = data.addElement("Gamma_Ydir_Err");
00044     i_gl_Gam_err    = data.addElement("Gamma_Err");
00045     
00046     i_gl_CsI_Xfit   = data.addElement("CsI_fitX");
00047     i_gl_CsI_Yfit   = data.addElement("CsI_fitY");
00048     i_gl_CsI_BD     = data.addElement("CsI_Bnd_Dst");
00049     i_gl_CsI_fitEr  = data.addElement("CsI_Fit_err");
00050     i_gl_CsI_fitErN = data.addElement("CsI_Fit_errNrm");
00051     i_gl_CsI_Xcount = data.addElement("CsI_Xtal_Count");
00052     i_gl_CsI_Xratio = data.addElement("CsI_Xtal_Ratio");
00053 
00054     
00055     compCode = 0;
00056 //THB    PrintControl::instance()->addPrinter("Recon comp code", new Printer_T<GlastRecon>(this)); 
00057     clear(); // just in case!
00058 
00059     
00060 
00061 }
00062 
00063 
00064 GlastRecon::~GlastRecon ()
00065 {
00066     clear(); // ensure all objects deleted
00067     delete mc;
00068     delete cal;
00069     delete trkr;
00070     delete trig;
00071     delete veto;
00072     delete prof;
00073 }
00074 
00075 void GlastRecon::accept (ReconVisitor& rv)
00076 {
00077     //  First do depend reconstructions
00078     Recon::accept(rv);
00079     
00080     //  Now do this one..
00081     rv.visit(this);
00082 }
00083 
00084 void GlastRecon::reconstruct (const GlastData* glastData)
00085 {
00086     if(state() == DONE) return;
00087     
00088     // no further analysis if no tracks 
00089     if(!trkr->xtrackList.size() || !trkr->ytrackList.size()) return;
00090 
00091     
00092     //   Compare Tracker Results to MC
00093     Vector t1 = Vector(trkr->data.getData("Fit_x_dir"),
00094                       trkr->data.getData("Fit_y_dir"),
00095                       trkr->data.getData("Fit_z_dir"));
00096     Vector t0 = Vector(trkr->data.getData("Gamma_x_dir"),
00097                       trkr->data.getData("Gamma_y_dir"),
00098                       trkr->data.getData("Gamma_z_dir"));
00099     Vector mct0 = Vector(mc->data.getData("MC_xDir"),
00100         mc->data.getData("MC_yDir"),
00101         mc->data.getData("MC_zDir"));
00102         if (trkr->data.getData("Fit_Type") == 0) {
00103                 t1=mct0;
00104                 t0=mct0;
00105         }
00106     data.addData(i_gl_xfit_err, t1.x() - mct0.x() );
00107     data.addData(i_gl_yfit_err, t1.y() - mct0.y() );
00108     data.addData(i_gl_zfit_err, t1.z() - mct0.z() );
00109     data.addData(i_gl_xGam_err, t0.x() - mct0.x() );
00110     data.addData(i_gl_yGam_err, t0.y() - mct0.y() );
00111     
00112     double mc_dot = t0*mct0;
00113     if(mc_dot > 1.0) mc_dot = 1.0; // protection for acos()
00114     else if(mc_dot <= -1.0 ) mc_dot = -0.999999;
00115     data.addData(i_gl_Gam_err, mc_dot<1? acos(mc_dot): 0);
00116     
00117     // now require energy (5 MeV?) in calorimeter
00118     if( cal->data["CsI_Energy_Sum"] < 0.005) return; 
00119 
00120     compCode = 1;
00121 
00122     // require track reconstruction to continue
00123     if (trkr->data["Fit_Type"] == 0) return;
00124 
00125     //   *** Some varibles to do analysis with *****
00126     float CsICorrEnergy = trkr->data["CsI_Corr_Energy"];
00127 
00128     // Tracker - calorimeter impact and matching varibles
00129     Point xCal = Point(cal->data.getData("CsI_X"),
00130         cal->data.getData("CsI_Y"),
00131         cal->data.getData("CsI_Z"));
00132     Point x0   = Point(trkr->data.getData("Fit_x0"), 
00133                       trkr->data.getData("Fit_y0"),
00134                       trkr->data.getData("Fit_z0"));
00135     
00136     float showerDist = (xCal.z()-x0.z())/t1.z();
00137     Point fitCsI  = x0 + t1*showerDist;
00138     data.addData(i_gl_CsI_Xfit, fitCsI.x() );
00139     data.addData(i_gl_CsI_Yfit, fitCsI.y() );
00140     
00141     float halfWidthX = .5*(DetectorData::Glast::xNum*(DetectorData::Tower::mod_width + .3) -.3);
00142     float halfWidthY = .5*(DetectorData::Glast::yNum*(DetectorData::Tower::mod_width + .3) -.3);   
00143     float CsIBndDistX = halfWidthX - fabs(fitCsI.x());
00144     float CsIBndDistY = halfWidthY - fabs(fitCsI.y());
00145     float CsIBndDist  = (CsIBndDistX < CsIBndDistY) ? CsIBndDistX : CsIBndDistY;
00146     data.addData(i_gl_CsI_BD, CsIBndDist );
00147     
00148     double d_par = t0*(x0 - xCal); 
00149     float CsIFitErr  = sqrt((x0-xCal).mag2() - d_par*d_par); 
00150     data.addData(i_gl_CsI_fitEr, CsIFitErr);
00151     
00152     //  Error model to normals impact error: a + b/E  ... why not a _+_ b/sqrt(E) ??? 
00153     float eFactor = .26 + .35 / CsICorrEnergy;
00154     float CsIFitErrNrm = CsIFitErr/eFactor;
00155     data.addData(i_gl_CsI_fitErN, CsIFitErrNrm );
00156     
00157     // Section to estimate number of hit CsI Xtals
00158     float num_layers = DetectorData::Tower::num_trays -1 - trkr->data.getData("fst_X_Lyr");
00159     num_layers /= DetectorData::Tower::num_trays -1.; 
00160     float xtal_count = 
00161         (CsICorrEnergy > 0)?10. + 2.7*log(CsICorrEnergy/.1) * (1. + .15*num_layers) : 0;
00162     float count_ratio = (xtal_count>0)? cal->data.getData("CsI_No_Xtals_Trunc")/xtal_count : 99.;
00163     data.addData(i_gl_CsI_Xratio, count_ratio);
00164     data.addData(i_gl_CsI_Xcount, xtal_count);
00165     
00166     // Section to set completion code
00167     setState(DONE);
00168     
00169     //      Minimal veto cut to trunc. cosmic runs
00170     if(    trig->data.getData("Trig_1x1") >  .5
00171         && trig->data.getData("Trig_1x1") < 5.5       
00172         && trig->data.getData("Trig_Level2_Bits") < 1.5 
00173         ) compCode = 2;
00174     
00175     //     Cuts to kill cosmics
00176     if(compCode == 2
00177         &&  trkr->data.getData("Surplus_Hit_Ratio") > 2.05
00178         &&  veto->data.getData("Veto_DOCA") > 30.
00179         &&        data.getData("CsI_Fit_errNrm") < 5.
00180         &&        data.getData("CsI_Xtal_Ratio") > .25
00181         ) compCode = 3;
00182     
00183     //    Cuts to improve resolution
00184     if(   compCode == 3 
00185         &&  trkr->data.getData("fst_Hit_Count") > 1.5
00186         &&  trkr->data.getData("fst_Hit_Count") < 2.8
00187         &&  trkr->data.getData("Diff_1st_XY_Lyr") > -.5
00188         &&  trkr->data.getData("Diff_1st_XY_Lyr") < .5
00189         &&  trkr->data.getData("Active_Dist") > 0.
00190         &&  trkr->data.getData("fst_X_Lyr") < 11.
00191         ) compCode = 4;
00192     
00193     //    Special Study HALT code 
00194     // if(   compCode == 4 
00195     //    &&  (trkr->data.getData("fst_Hit_Count") <= -1 ||
00196         //      trkr->data.getData("fst_Hit_Count") >= 4)
00197     //    ) compCode = 5;
00198     if(compCode >=4 && data.getData("Gamma_Err") > 0.12)
00199                 compCode = 5;
00200 
00201 
00202 }
00203 
00204 void GlastRecon::clear ()
00205 {
00206     Recon::clear();
00207     data.clear();
00208     compCode = 0;
00209 }
00210 void GlastRecon::printOn(std::ostream& out)const
00211 { 
00212     out << "Completion code:" << completionCode() << std::endl;
00213 }
00214 
00215 const char* GlastRecon::nameOf () const{    return "GlastRecon"; }